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Decomposed  Analysis  of  the  Reliability  of  Area-Array  Electronic  Packages 
Including  Assembly  Process  Effects 

Thesis  directed  by  Associate  Professor  Ganesh  Subbarayan 

The  research  presented  in  this  thesis  focuses  on  the  use  of  a  decomposed 
analysis  technique  to  assess  the  reliability  impact  on  area-array  packages  due  to 
assembly  process  effects.  The  developed  models  enable  one  to  determine  the 
equilibrium  configuration  of  electronic  packages  with  warped  and  non-warped  solder 
joints  with  significant  computational  efficiency  at  a  reasonable  accuracy.  The  code 
described  herein  uses  a  nonlinear  optimization  procedure  that  ensures  the 
approximate  satisfaction  of  the  energy  balance  equation.  The  developed  procedure 
is  demonstrated  on  two-  and  three-dimensional  hypothetical  and  “real-world” 
electronic  packages  with  warped  and  non-warped  solder  joints.  It  is  shown  that  with 
the  use  of  the  decomposed  solution  methodology,  for  a  225  I/O  PBGA  package  with 
non-warped  solder  joints,  a  speedup  of  nearly  seven  times  is  achieved  at  an  accuracy 
loss  in  displacements  of  approximately  5.46%.  Since  the  analysis  procedure  is 
independent  of  the  number  of  solder  interconnects,  significantly  larger  timesavings 
are  expected  for  larger  packages. 

One  ultimate  utility  of  any  technique  employed  to  analyze  electronic 
packages  is  to  assess  package  reliability.  Reliability  assessment  techniques  normally 
consist  of  three  basic  steps  to  calculate  (Subbarayan,  et  al.,  1996): 

1.  solder  joint  shapes. 


IV 


2.  thermal  stress/strain  distribution  in  the  most  susceptible  joint  after 
solidification  and  during  thermal  cycling,  and 

3.  fatigue  lives  of  the  most  susceptible  or  any  other  selected  solder 
joints. 

For  step  one,  we  will  demonstrate  a  novel  methodology  for  studying  the 
effect  of  printed  circuit  board  (PCB)  warpage  and  solder  joint  volume  variation  on 
the  final  equilibrium  configuration  of  area-array  packages.  The  effect  of  warpage  is 
analyzed  using  a  two-step  procedure.  In  the  first  step,  the  restoring  forces  and 
moments  (in  the  molten  state  of  solder  droplet)  that  result  from  a  given  solder  joint 
height,  solder  material  volume,  pad  diameter,  and  pad  tilt  are  predicted  using  the 
surface  tension  theory.  In  the  second  step  of  the  analysis,  the  forces  and  moments  at 
individual  solder  joints  caused  by  varying  solder  heights,  pad  tilts,  and  solder 

volume  are  combined  using  an  optimization  procedure  to  predict  the  equilibrium 
configuration  of  the  package. 

The  developed  procedure  is  demonstrated  on  a  hypothetical  area-array 
package  with  nine  solder  joints.  Through  an  analysis  of  two  scenarios  for  the 
hypothetical  array:  1)  constant  solder  volume  with  symmetric  and  non-symmetric 
package  warpage,  and  2)  linearly  distributed  solder  volume  without  warpage,  it  is 
shown  that  printed  circuit  board  warpage  can  cause  packages  to  tilt  during  solder 
reflow  resulting  in  variations  in  solder  joint  heights  and  pad  tilts. 

For  step  two  in  assessing  package  reliability,  namely  thermal  stress/strain 
distribution,  we  will  combine  the  decomposition  based  analysis  technique  and  the 
methodology  for  predicting  the  equilibrium  configuration  of  warped  electronic 
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packages  to  study  the  reliability  impact  of  warped  solder  joints  in  area-array 
packages.  The  methodology  is  developed  and  demonstrated  on  a  hypothetical  area- 
array  package  with  four  warped  solder  joints.  An  assumption  that  enabled 
computational  efficiency  was  that  the  geometrical  change  (to  the  nominal  geometry) 
caused  by  warpage  of  the  component/PCB  did  not  significantly  alter  the  stiffness  of 
the  structure.  This  allowed  the  use  of  substructure  models  of  the  component/PCB 
together  with  accurate  models  of  individual  solder  joints  with  widely  varying  shapes 
in  the  decomposed  analysis  procedure. 

It  is  shown  that  the  use  of  the  decomposed  solution  methodology  for 
analyzing  a  package  with  warped  solder  joints  leads  to  an  overall  accuracy  loss  in 
relative  shear  displacements  of  approximately  7. 18  percent.  This  error  translates  to 
an  accuracy  loss  in  fatigue  life  calculation  of  approximately  1 1.83  percent.  A 
comparison  of  the  reliability  of  a  package  with  warped  and  non-warped  joints 
reveals  package  tilt  due  to  solder  volume  variation  and/or  non-symmetric  warpage 
can  change  the  stand-off  height  of  solder  joints  profoundly  affecting  solder  joint 
reliability.  However,  we  found  that  individual  solder  pad  tilts  in  a  range  from  0  to 
3  degrees,  independent  of  solder  standoff  height,  do  not  affect  the  reliability  of 
solder  joints. 
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Chapter  1 


Introduction 


1.1  Decomposition 

Engineering  systems,  in  general,  are  physically  challenging  to  model  and 
computationally  expensive  to  simulate.  Electronic  packages  are  excellent  examples 
of  complex  engineering  systems  since  they  involve  multiple  physical  domains, 
namely  electrical,  mechanical,  and  thermal  and  complexity  within  each  domain.  The 
analysis  complexity  arises  out  of  factors  such  as:  nonlinear  behavior  of  materials 
(e.g.,  plastic  and  creep  deformations  of  solder),  complex  package  construction,  and 
the  large  number  of  interconnections.  As  electronic  packages  increase  in  size  and 
density,  common  analysis  techniques  such  as  finite  element  (FE)  analysis  have 
become  very  inefficient  for  these  systems  and  are  too  slow  for  quick  design 
decisions.  Optimal  design  of  electronic  packages  is  even  more  expensive  since 
optimization  techniques  are  iterative  in  nature  and  each  iteration,  using  FE  analysis, 
would  require  regeneration  of  the  FE  mesh  (for  the  new  geometry)  prior  to  re¬ 
analysis.  The  need  for  efficient  as  well  as  accurate  techniques  for  analysis  and 
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design  of  complex  mechanical  systems  in  general  and  electronic  packages  in 
particular  are  obvious. 

A  brief  background  on  the  methodology  developed  by  Deshpande  and 
Subbarayan  (1999)  as  it  relates  to  the  present  research  is  described  below.  In 
general,  the  finite  element  method  is  a  popular  choice  for  the  analysis  of 
microelectronic  packages.  The  concept  of  a  superelement  is  well  established  in  the 
finite  element  method,  where  a  global  stiffness  matrix  is  partitioned  into  retained  and 
eliminated  portions  (referred  by  subscripts  r  and  e  below)  as  follows: 


K 

K 


rr 

re 


KJk 

kJK 


1.1 


K„ur=fr 

Krr=Krr-KreKe;'Ker  1.2 

fr=fr-KreKeelfe 

The  concept  of  superelements  in  finite  element  analysis  allows  the 
partitioning  of  a  large  linear  structure  into  substructures  for  ease  of  analysis  (see  for 
example,  Zienkiewicz  and  Taylor,  1990).  When  two  superelements  are  “connected”, 
the  mesh  or  the  nodal  degrees  of  freedom  at  the  interface  must  match.  Also, 
superelements  are  valid  only  for  linear  systems  since  they  rely  on  the  idea  of  a 
constant  stiffness  matrix.  Thus,  their  applicability  to  microelectronic  packages  is 
limited  since  solder  alloys  and  other  materials  commonly  used  are  known  to  behave 
nonlinearly.  The  decomposition  methodology  extends  the  concept  of  superelement 
by  allowing  for  mesh  mismatch  between  the  substructures  at  the  interface,  and  by 
allowing  nonlinearly  behaving  substructures. 
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The  earlier  derived  theoretical  methodology  of  Deshpande  and  Subbarayan 
(1999)  were  validated  in  their  work  using  a  code  developed  specifically  for  the 
analysis  of  the  5x5  array  used  in  that  study.  Any  change  in  package  configuration  or 
construction  requires  a  considerable  amount  of  reprogramming  and  debugging.  The 
application  and  demonstration  of  a  code,  based  on  the  decomposition  procedure  and 
assumed  displacements  fields,  to  hypothetical  and  “real-world”  electronic  packages 
is  one  of  the  goals  of  the  current  research.  By  simply  modifying  the  code  input  deck 
to  account  for  new  package  geometries,  we  were  able  to  analyze  new  electronic 
packages  with  ease. 

1.2  Assembly  Process-Induced  Variations  in  Array-Packages 

The  current  trends  in  the  electronic  packaging  industry  are  toward  larger 
packages,  increased  density,  new  substrate  materials  and  new  printed  circuit  board 
(PWBs)  constructions.  This  has  led  to  issues  with  PWB  and  component  warpage 
which  is  believed  to  lead  to  component  misalignment  and/or  tilt  during  solder  reflow 
processes  resulting  in  a  premature  solder  joint  failure  during  use  conditions.  There  is 
a  need  to  systematically  study  the  reliability  implications  of  symmetric  and  non- 
symmetric  circuit  board/component  warpage  and  solder  volume  variation  in 
electronic  packages  since  it  is  believed  that  these  variations  could  lead  to  non- 
intuitive  failures  in  joints  other  than  the  usual  “corner  joint”  at  the  die  edge. 

1.3  Literature  Review 


1.3.1  Decomposition  in  Electronic  Packages 
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Deshpande  accomplished  much  of  the  literature  review  on  domain  and  design 
composition  during  his  Ph.  D.  research.  The  application  of  decomposition-based 
finite  element  analysis  to  electronic  packages  is  novel  work  and  Deshpande  (1999) 
presents  the  only  literature  found  on  the  application  of  these  techniques  for  efficient 
analysis  of  electronic  packages. 

Cheng,  et  al.,  (1998)  present  an  efficient  method  for  FE  analysis  of  BGA 
packages.  Their  method  is  based  on  replacing  solder  joints  by  equivalent  beams  and 
calculating  the  properties  of  the  equivalent  beam  using  a  designed  experiment  based 
on  analyses  of  3-D  solder  models.  Once  the  equivalent  beam  is  formulated  they 
perform  a  full  3-D  FE  analysis  of  the  entire  package  using  the  equivalent  beam.  The 
methodology  is  not  decomposition  based  since  the  component  and  PWB  are  not 
separate  subsystems,  thus  any  changes  in  either  would  require  a  complete  re-analysis 
of  the  3-D  model  with  the  beams.  The  decomposition  methodology  is  free  from  this 
limitation. 

Saito,  et  al.,  (1999)  also  present  an  efficient  stress  and  displacement  analysis 
method  for  area-array  structures  based  on  the  equivalent  beam  methodology  where 
two-node  beam  elements  with  stiffness  matrices  extracted  from  load-displacement 
relationships  are  used  to  model  the  solder  joints  and  shell  elements  are  used  to  model 
the  package  and  PWB.  This  method  has  the  same  limitation  as  Cheng’s  method 
above. 

Darbha  and  Dasgupta  (1999)  present  a  stress  analysis  technique  for  electronic 
packages  based  on  a  nested  finite  element  methodology  (NFEM)  where  colonies  of 
nested  sub-elements  are  created  inside  a  main  element  to  capture  localized  large 
gradients  of  stress  and  strain.  The  nesting  concept  reduces  the  number  of  degrees  of 
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freedom  of  the  FE  model,  thus  reducing  computational  time  since  the  nested  finite 
elements  are  only  developed  for  critical  regions  of  the  package.  Once  again,  this 
methodology  is  not  decomposition  based  since  any  design  change  would  require  re¬ 
analysis  of  a  full  3-D  finite  element  mesh. 

1,3.2  Warpage  in  Electronic  Packages 

A  large  body  of  literature  is  published  on  warpage  in  electronic  packages 
and  its  effect  on  solder  joint  reliability.  Most  of  the  literature  is  aimed  at  doing  one 
of  three  things: 

1.  Development  of  methodologies  to  predict  warpage  of  electronic  packages 
with  or  without  analysis  to  determine  solder  joint  reliability. 

2.  Development  of  techniques  to  provide  “real-time”  monitoring  of  warped 
electronic  packages  pre,  in-situ,  or  post-assembly  in  order  to  determine 
solder  yield,  where  the  yield  criteria  is  based  on  a  maximum/minimum 
solder  joint  height  that  would  cause  an  open  or  short  circuit. 

3.  Development  of  techniques  to  predict  the  shape  and  yield/reliability  of 
solder  joints  in  warped  electronic  packages. 

Amagai  (1999)  found  that  package  warpage  tremendously  affected  solder 
joint  reliability  in  Chip  Scale  Packages.  The  effect  of  molding  compound  and  die 
attach  film  on  package  warpage  was  investigated  using  3-D  non-linear  (viscoplastic 
and  viscoelastic)  finite  element  analysis.  Amagai  was  able  to  correlate  the 
magnitude  of  package  warpage  after  encapsulation  to  solder  fatigue  life.  Thermal 
cycling  of  the  3-D  finite  element  models  and  subsequent  fatigue  life  calculations 
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showed  the  solder  joints  near  the  chip  edge  of  the  package  with  the  greatest  warpage 
failed  first. 

Yao  and  Qu  (1999)  investigated  the  effect  of  printed  wiring  board  size  on  the 
overall  deflection  of  flip-chip  assemblies  and  found  that  the  size  of  the  PWB  had 
significant  implications  on  the  design  of  multi-chip  modules  (MCMs).  Three- 
dimensional  finite  element  analyses  were  used  to  study  the  stress  distribution  in  and 
deflection  of  the  flip  chip  assembly.  The  FE  model  was  validated  with  experimental 
results  of  stress  distribution.  The  results  indicate  that  a  square  array  placement 
pattern  is  preferable  to  a  staggered  array  for  MCMs  and  a  minimum  distance 
between  adjacent  dies  should  be  approximately  2  times  the  die  size  in  order  to 
reduce  mechanical  interaction  between  chips. 

He,  X.  and  Lu,  S.  (1998)  present  and  demonstrate  a  new  portable  projected 
grating  moire  system.  The  system  was  designed  and  built  for  component  level 
warpage  measurement  and  allows  precision  measurement  of  warpage,  providing 
valuable  insight  into  package  reliability,  manufacturing  process  and  model 
validation.  Verma,  et  al.,  (1998)  developed  a  compact  apparatus  combining  far 
infra-red  Fizeau  interferometry  and  shadow  moire  with  enhanced  sensitivity 
methodologies  to  measure  real-time  warpage  with  a  variable  sensitivity  ranging  from 
2.5  |im  to  100  pm. 

Lee  (1997)  developed  an  analytical  model  to  the  study  warpage  of  plastic 
ball-grid-array  packages  under  uniform  thermal  loading.  The  model  was  derived 
from  the  classical  laminated  plate  theory.  A  series  of  parametric  studies  were 
performed  to  investigate  the  effects  of  chip  size,  chip  thickness,  BT  substrate 
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thickness,  and  molding  compound  thickness  on  package  warpage.  The  model 
provides  the  packaging  engineer  with  a  convenient  tool  to  evaluate  PBGA  warpage. 

Egan,  et  al.,  (1998)  described  the  overall  warpage  of  a  PBGA  package  as 
consisting  of  two  curvatures,  a  die-area  curvature  and  a  mold-area  curvature.  This 
was  based  on  measurements  done  on  physical  samples  of  PBGAs.  An  analytical 
expression  for  dual-curvature  beam  bending  is  developed  and  compared  with 
previous  predictions  for  multi-layered  beam  bending  and  3-D  finite  element  analysis. 
The  dual-curvature  approach  yielded  very  good  2-D  accuracy  and  fairly  good  3-D 
accuracy. 

In  all  of  the  above  research,  solder  joint  shape  prediction  and  the  reliability 
resulting  from  warped  circuit  boards  and  packages  has  not  received  much  attention. 
Ju  et  al.,  (1994),  Chan  et  al.,  (1997),  and  Tower,  et  al.,  (1999)  appear  to  be  the  only 
studies  of  this  nature. 

Ju  et  al.,  (1994)  studied  the  effects  of  manufacturing  variations  on  the 
reliability  of  solder  joints  between  a  ceramic  ball  grid  array  (CBGA)  electronic 
package  and  PWB.  Manufacturing  variations  considered  included  changes  in  solder 
volume,  joint  height,  and  pad  size,  where  the  joint  height  variation  was  caused  by  the 
package  warpage.  They  investigated  the  effect  of  spacers  to  model  warpage  and 
control  solder  joint  height  and  compared  the  fatigue  life  of  the  package  outermost 
comer  joint’s  convex  (truncated  sphere  shaped),  cylindrical  or  concave  (hour-glass 
shaped)  profiles  produced  by  these  spacers.  Solder  life  prediction  calculations 
showed  concave  solder  joints  had  long  fatigue  lives  and  were  less  sensitive  to 
manufacturing  variations  when  the  spacers  were  used.  They  concluded  the 
reliability  of  convex  solder  joints  is  significantly  affected  by  the  height  variations 
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associated  with  the  package  warpage.  More  than  50  percent  variation  in  lives 
estimated  for  a  +/-  10  percent  height  change. 

Chan,  et  al.,  (1997)  established  a  reliability  model  to  study  thermal  fatigue  in 
PBGA  assemblies.  The  model  was  used  to  simulate  a  configuration  with  a  large 
number  of  warpage  affected  solder  joints.  A  regression  model  based  on  solder 
restoring/reaction  forces  during  reflow  was  used  to  predict  the  equilibrium 
configuration  of  a  warped  electronic  package.  Analyzing  two  PBGA  electronic 
packages,  a  72  I/O  cavity-up  and  a  540  I/O  cavity-down  assembly,  the  study  showed 
the  reliability  of  the  cavity-down  assembly  was  strongly  affected  by  warpage.  Arch 
and  saddle-type  warpage  changed  the  shape  of  the  comer  solder  joint  by  reducing  its 
height,  and  thus,  fatigue  life  considerably.  The  study  inferred  a  bowl-shaped 
warpage  would  increase  the  life  of  the  comer  joint  beyond  that  of  the  no-warpage 
model  because  the  joint’ s  height  would  be  increased. 

Tower,  et  al.,  (1999)  developed  a  model  to  predict  yield  of  flip-chip  solder 
assemblies  during  solder  reflow  by  considering  mean  and  standard  deviation  of 
volume  distribution,  the  assembly  warpage,  the  die  size,  and  the  number  of 
input/outputs.  Again,  a  regression  model  based  on  solder  restoring/reaction  forces 
during  reflow  was  used  to  predict  the  equilibrium  configuration  of  a  warped 
electronic  package.  The  model  showed  a  strong  correlation  between  warpage  and 
assembly  yield,  as  warpage  increased,  yield  decreased. 

The  above-cited  studies  have  ignored  the  pad  rotation  that  would  naturally 
result  from  warped  electronic  packages.  Renken  and  Subbarayan  (2000)  have  shown 
that  a  5°  pad  rotation  could  nearly  halve  the  standoff  height  in  area-array  solder 
joints  due  to  the  resulting  nonsymmetric  distribution  of  droplet  volume.  Thus,  the 
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inclusion  of  pad  rotation  in  the  present  study  is  believed  to  have  important 
consequences  for  the  solder  joint  reliability  and  the  belief  that  warpage  could  lead  to 
non-intuitive  failures  in  joints  other  than  the  usual  “comer  joint”  at  the  die  edge. 

1.4  Summary  of  Contributions  and  Thesis  Outline 

1.4.1  Contributions 

The  major  contributions  of  thesis  can  be  broken  into  three  areas: 

1.  An  implementation  procedure  for  the  analysis  of  electronic  packages 
based  on  the  decomposition  methodology  enabling  one  to  determine  the 
equilibrium  configuration  of  electronic  packages  with  significant 
computational  efficiency  at  a  reasonable  accuracy. 

2.  Development  of  a  model  and  an  accompanying  computer  code  for 
studying  the  effect  of  printed  circuit  board  (PCB)/package  warpage  and 
solder  volume  variation  on  the  final  equilibrium  configuration  of  area- 
array  packages  during  solder  reflow. 

4.  Development  of  a  methodology  for  the  analysis  of  electronic  packages 
with  warped  solder  joints  based  on  the  decomposition  procedure  enabling 
one  to  carryout  a  systematic  study  of  the  reliability  implications  of 
warpage  in  electronic  packages. 

1.4.2  Thesis  Outline 

In  Chapter  2  we  introduce  the  decomposition  procedure  for  the  analysis  of 
electronic  packages  and  the  code  developed  utilizing  this  procedure  to  analyze 
electronic  packages.  The  chapter  concludes  with  the  developed  procedure 
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demonstrated  on  two-  and  three-dimensional  hypothetical  and  “real-world” 
electronic  packages  with  non-warped  solder  joints. 

In  Chapter  3  we  describe  a  model  and  a  computer  code  for  studying  the  effect 
of  printed  circuit  board  (PCB)/package  warpage  and  solder  volume  variation  on  the 
final  equilibrium  configuration  of  area-array  packages  during  solder  reflow  process. 
The  chapter  concludes  with  the  developed  procedure  demonstrated  on  hypothetical 
electronic  packages  with  warped  solder  joints. 

In  Chapter  4  we  demonstrate  the  decomposition  procedure  for  analyzing 
electronic  packages  with  warped  solder  joints.  We  then  demonstrate  the  procedure 
on  a  three-dimensional  hypothetical  area-array  electronic  package.  We  conclude  the 
chapter  with  a  discussion  on  the  reliability  impact  of  electronic  packages  with 
warped  solder  joints  and  the  consequences  of  pad  rotation/tilt  on  solder  joint 
reliability. 

Chapter  5  concludes  the  thesis  by  highlighting  the  novel  contributions  in  the 
present  study  and  by  making  recommendations  for  future  work  for  further 
continuation  of  the  present  study. 


Chapter  2 


Applications  of  a  Decomposed  Analysis  Procedure  for  Area-Array 

Packages 


2.1  Chapter  Overview 

In  this  chapter  we  introduce  the  decomposition  procedure  for  analysis 
of  electronic  packages  and  the  code  developed  utilizing  this  procedure  to 
analyze  electronic  packages.  The  developed  procedure  is  demonstrated  on 
two-  and  three-dimensional  hypothetical  and  “real-world”  electronic 
packages  with  non-warped  solder  joints. 

Finally  to  illustrate  the  usefulness,  flexibility  and  strength  of  this 
procedure,  we  review  displacement  accuracies,  computational  efficiencies, 
and  superelement  energy  contributions. 

The  objectives  of  this  chapter: 

•  Refine  model  for  solder  response  through  additional  basis  terms. 

•  Demonstrate  methodology  on  “real-world”  electronic  package. 

•  Provide  insight  into  decomposed  approach. 

•  Displacement  versus  relative  displacement  errors. 


26 


•  Provide  insight  into  superelement  energy  contributions  to  relative 
displacement  errors. 

2.2  Introduction 

New  products  are  pushing  density  and  speed  of  electronic  packages  to  the 
limits.  To  support  this  evolution,  packaging  technologies  are  moving  at  a  rapid  pace. 
It  is  believed  by  the  year  2012,  complex  electronic  packages  will  contain  over  5000+ 
I/Os  (SIA,  1997).  Analysis  techniques  must  evolve  at  an  even  faster  pace  if  quick 
decisions  regarding  package  reliability  are  to  be  made.  Numerical  techniques  such 
as  finite  element  analysis  are  commonly  used  in  determination  of  stresses  and 
strains,  and  through  these  the  package  reliability.  However,  with  decreasing  package 
size,  increasing  I/O  count,  and  complex  package  configurations  and  materials,  finite 
element  analysis  has  become  computationally  expensive.  Analysis  of  packages  can 
take  hours  or  even  days  to  complete  at  the  present  time,  owing  primarily  to  the 
nonlinear  behavior  of  materials. 

Deshpande  and  Subbarayan  (1999)  developed  and  demonstrated  an  efficient 
technique  for  the  design/analysis  of  electronic  packages.  This  procedure  is  inspired 
by  the  domain  and  design  decomposition  methodologies  literature.  The  former 
methods  largely  address  the  issue  of  partitioning  a  large  domain  (approximated  by  a 
finite  element  mesh)  into  well-balanced  subdomains  (submeshes)  that  can  be 
efficiently  solved  on  a  parallel  computer  (see  for  example  Farhat  et  al.  1987,  1989). 
The  latter  methods  aim  to  divide  a  large  design  problem  (often  posed  as  an 
optimization  problem)  into  a  series  of  simpler  subproblems  for  solution  efficiency 
(Haftka  and  Gurdal,  1992).  The  procedure  of  Deshpande  and  Subbarayan  (1999)  is 
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among  the  first  for  simultaneous  design  and  domain  decomposition.  They 
demonstrated  a  speed-up  of  approximately  350%  at  an  accuracy  loss  of  only  6%  on  a 
hypothetical  5x5  area  array  procedure.  Their  analysis  of  the  floating  point 
operations  demonstrated  an  unbounded  increase  in  speed-up  for  larger  array  sizes. 

During  the  present  research,  we  develop  a  general  implementation  scheme 
for  the  technique  developed  by  Deshpande  and  Subbarayan.  The  developed  code 
addresses  a  critical  need  in  the  electronic  industry  for  quick  reliability  decisions. 
One  major  aspect  of  the  proposed  implementation  scheme  is  the  parametric 
representation  of  the  overall  problem  that  allows  the  system  response  to  be 
determined  by  a  relatively  few  unknowns.  This  is  in  contrast  to  the  nonlinear  finite 
element  analysis  schemes  that  would  require  an  iterative  solution  to  a  large  number 
of  nodal  unknowns.  A  second  major  aspect  of  the  implementation  scheme  is  the 
underlying  optimization  procedure  that  iteratively  determines  the  solution  to  the 
unknowns  for  any  area  array  package  described  in  the  required  manner. 

2.3  Methodology 

The  first  step  in  the  present  application  of  the  procedure  is  a  partitioning  of 
the  package  into  solder  joints,  the  component,  and  the  circuit  board  (Figure  2.1). 
Further,  the  parameters  that  describe  the  package  geometry  and  its  deformation  are 
categorized  as  design  variables  and  analysis  variables.  The  variables  that  determine 
the  geometry  of  the  package,  or  the  design  variables,  are:  the  sizes  of  the  package 
and  the  PCB,  the  pad  diameters,  the  solder  volume  and  the  weight  of  the  package 
borne  by  solder  joint.  The  analysis  variables  are  typically  the  displacements 
determined  for  a  given  geometry.  These  displacements  could  occur  either  within 
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each  subsystem  (solder  joint,  PCB,  and  component)  or  at  the  interface  between  the 
solder  joints  and  the  component  as  well  as  between  the  solder  joints  and  the  PCB. 
These  latter  displacements  occurring  at  the  interface  are  in  some  sense  more  critical 
since  these  displacements  characterize  the  interaction  between  the  subsystems  joined 
at  the  interface  (Deshpande  and  Subbarayan,  1999). 


PCB  Package/  Solder 

Component 

Figure  2.1  Decomposition  of  electronic  packages. 


2.3.1  Problem  Parameterization 

The  novelty  of  the  proposed  implementation  scheme  lies  in  the  assumption  of 
displacement  fields  at  solder/component  and  solder/PCB  interfaces  and  the  ability  of 
the  code  to  generate  these  displacements  for  any  type  of  array  package  based  on  a 
general  code  input  deck.  Each  individual  solder  joint  is  subjected  to  relative  shear 
(u-v)  and  axial  (w)  displacements.  These  displacements  are  due  to  the  global 
mismatch  in  the  values  of  the  coefficients  of  thermal  expansion  for  the  PCB  and  the 
component.  In  addition,  the  solder  joint  is  subjected  to  u,  v,  and  w  displacements 
due  to  the  local  mismatch  in  thermal  expansions  between  the  solder  and  either  the 
PCB  or  the  component.  In  general,  the  u  ,v,  and  w  displacements  across  a  solder 
joint  due  to  the  combination  of  global  and  local  mismatches  may  be  non-uniform. 
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As  a  first  order  approximation,  we  assume  this  variation  to  be  linear.  The 
differential  expansion  due  to  the  local  mismatch  at  any  point  on  the  interface 
between  the  solder  joint  and  the  component  or  PCB  is  assumed  radial  in  direction 
with  respect  to  the  center  of  the  solder  joint.  Based  on  our  assumption  of  linear 
variations  in  local  displacements,  the  displacement  field  across  the  top  and  bottom 
solder  joint  interfaces  can  be  approximated  by  eight  variables,  £o,  ,  £2,  and  £3  at  the 

bottom  interface,  and  yo.  Yi  ,  Y2,  and  73  at  the  top  interface  as  shown  in  Figure  2.2. 


Node  Analyzed 


PCB/Solder  Interface 

=>  Ux  ^  *  a,  K  B  *  AT  *  Rnp  *  CosOjO  + 

5,  *  AbsfO*  K'R  -  ocx  So,'fcr)  *  AT  *  Rnodc  *  Cos(0) 

=>  Uy  ^  ^  ^  *  AT  *  Rnp  *  Sin(v)  + 

*  (Oy  R  Oy  ***)  *  AT  *  Rmide  *  Sin(6) 

=>  U,wlJc#  =  *  a,  tK  B  *  AT  +  *  Rnodc  *  Cos(<|>) 

Component/Solder  Interface 
=*  U„ =  ft,  *  O, PrB  *  AT  *  R„r  -  T„)  COS(V)  + 

Y,  *  Absfa,  r""r"“  -  a,  *  AT  *  *  Cos(e> 

=»  U,  “**  =  <&,  *  a, K'B  *  AT  *  R„r  -  Yo)  SinOlO  + 

Y  *  Absfa,  O, Sol,kr)  *  AT  *  *  Sin(6) 

=» U,“"M  =  \2*  a, K" t1™ *  AT  +  Y:  +  y,  *  *  Cos(0) 

5o,5„  Y2-Y3 are  unknown  parameters  at  solder  joint 


Figure  2.2  Generating  displacement  fields  in  the  package  decomposition. 


At  the  bottom  interface,  the  variable  £,0  accounts  for  the  u-v  displacements 
due  to  thermal  expansion  of  the  PCB,  whereas  the  variable  accounts  for  the  u-v 
displacements  due  to  local  mismatch  at  the  PCB/solder  interface  .  The  actual  nodal 


30 


u-v  displacements  are  calculated  by  adding  the  local  and  global  displacements.  The 
axial  displacements  are  assumed  to  occur  along  the  line  connecting  the  package 
neutral  point  to  the  solder  center  joint  and  is  represented  by  a  uniform  w 
displacement  accounted  for  by  variable  ^  and  a  linear  variation  of  slope 
superimposed  on  w  and  accounted  for  by  variable  The  actual  nodal  w 
displacements  are  calculated  by  adding  the  uniform  and  the  linearly  varying 
components  of  the  w  displacements.  Variable  £2  can  also  account  for  any 
bending/warpage  of  the  package  due  to  thermal  mismatch. 

At  the  top  interface,  the  variable  yo  accounts  for  the  relative  shear 
displacement  of  the  solder  joint  due  to  global  thermal  mismatch.  From  this,  the 
actual  nodal  u-v  displacements  of  the  component  are  calculated  by  subtracting  this 
relative  shear  displacement  from  the  PCB  u-v  displacement  and  adding  the  local  u-v 
interfacial  displacements,  accounted  for  by  variable  yi.  The  relative  axial 
displacement  of  the  solder  joint,  y2,  and  linear  varying  slope  imposed  on  this 
displacement,  73,  are  added  to  the  axial  displacements  of  the  PCB  to  generate  the 
axial  displacements  of  the  component  (see  Figure  2.2). 

The  use  of  the  package  temperature  rise  (AT),  thermal  expansion  coefficients 
(a),  and  geometric  entities,  such  as  PCB  thickness  (tPCB),  in  calculating  the  nodal 
displacements  serve  to  normalize  our  displacement  field  variables  making  the 
optimization  variables  scale  independent.  The  approximation  of  the  package 
behavior  using  these  eight  variables  greatly  reduces  the  total  number  of  unknowns  at 
each  interface,  compared  to  a  full  nonlinear  finite  element  analysis. 
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Deshpande  and  Subbarayan  (1999)  used  only  three  variables  per  solder 
interface  (six  variables  per  solder  joint)  to  approximate  the  relative  displacement 
field,  u,  v  and  6 ,  across  a  solder  joint.  The  displacement  u  is  the  relative  shear 
displacement  due  to  global  mismatch.  The  axial  displacements  are  represented  by  a 
uniform  displacement  v  and  a  linear  variation  with  slope  6  superimposed  on  v  (see 
Figure  2.3).  In  their  analysis,  the  local  displacements  are  assumed  constant  at  each 
interface. 


Figure  2.3  Relative  displacements  in  the  package  decomposition. 

In  our  analysis,  the  displacements  at  every  interfacial  node  are  captured  using 
four  unknown  variables,  for  a  total  of  eight  variables  per  solder  joint.  These 
variables  are  selected  based  on  physical  intuition  and  although  there  is  an  increase  in 
the  number  of  optimization  variables  per  solder  joint  (eight  versus  six),  the  present 
procedure  provides  a  more  structured  and  robust  approach  to  analyzing  electronic 
packages.  This  approximation  is  one  key  factor  in  increasing  the  efficiency  of  the 
decomposed  analysis  approach.  For  a  given  displacement  field,  the  response  of  the 
solder  joint  was  described  by  the  rate  at  which  work  is  done  by  the  solder  joint  on 
the  rest  of  the  system  and  vice  versa. 
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2.3.2  Subsystem  Interfacial  Work 

In  area-array  packages,  the  materials  used  for  the  construction  of  PCB  and 
the  component  can  often  be  approximated  as  linearly  behaving.  It  should  be  noted 
that  in  a  general  application  of  the  partitioning  methodology,  the  linear  substructure 
need  not  be  a  clearly  identifiable  subsystem  (such  as  the  PCB),  and  could  well  be  a 
linearly  behaving  subdomain  of  one  homogeneous  structure.  The  rate  at  which 
energy  is  stored  in  these  linear  substructures  (assumed  to  be  discretized  by  a  finite 
element  mesh)  for  any  arbitrary  displacement  u  at  the  retained  nodes  is  given  by: 

W  =  J  <7pq€pqdV  =  uTKu  2.1 

where  Gpq  is  the  stress  tensor,  e  is  the  rate  of  strain  tensor,  and  K  is  the  reduced 

stiffness  matrix  of  the  substructure  or  “superelement.”  The  superelements  were  built 
using  the  ABAQUS  finite  element  code  in  the  present  study,  but  most  finite  element 
codes  allow  the  calculation  of  these  matrices. 

Typically,  in  any  electronic  package,  the  solder  joints  are  nonlinearly 
behaving  subsystems  and  a  detailed  nonlinear  finite  element  analysis  is  necessary  to 
capture  the  response  of  these  joints.  In  the  developed  procedure,  the  interaction 
between  the  solder  joints  and  the  rest  of  the  system  (PCB  or  component)  at  an  instant 
t  is  captured  through  the  rate  at  which  work  is  done  at  the  solder  interface,  which  is 
calculated  as: 

W=\Tiuids^Fi{t)iii{t) 

S  '=! 


2.2 
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where  n  is  the  total  number  of  nodes;  Ft  and  ui  are  the  forces  and  velocities 
respectively,  at  node  i. 

2.3.3  Coordinating  Subsystems  Towards  System-level  Solution 

For  the  overall  energy  balance  of  the  system,  the  three  subsystems  mentioned 

earlier  must  each  satisfy  the  equations  governing  balance  of  mechanical  energy,  as 
must  the  package  as  a  whole.  Ignoring  the  rate  of  kinetic  energy  enables  one  to 
define  the  equilibrium  configuration  of  the  system.  The  error  introduced  due  to  the 
partitioning  of  the  package  can  be  broadly  divided  into  two  categories:  those  that 
accumulate  at  the  interface  between  two  subsystems  due  to  incorrectly  imposed 
displacements/tractions  at  the  interface,  and  those  that  occur  during  the  solution  of 
the  energy  balance  equation  by  the  finite  element  method.  Since  the  latter  error  is 
unavoidable  even  when  an  analysis  is  carried  out  without  any  decomposition,  the 
former  error  occurring  at  interfaces  is  the  one  that  is  caused  by  the  decomposition 
process.  Deshpande  and  Subbarayan  (1999)  derived  a  relationship  for  this  error  for  a 
general  nonlinear  system  loaded  over  a  step  lasting  a  time  At  as: 

K 

£  =  I  f  (W,+W,)<*  2.3 

k  t 

where  W  is  the  rate  of  external  work,  subscripts  p  and  q  refer  to  the  two  subsystems 
that  share  the  interface  k,  and  K  is  the  total  number  of  new  interfaces  introduced  due 
to  the  partitioning  process.  K  will  in  general  equal  one  less  than  the  total  number  of 
subsystems.  If  the  tractions  and  displacements  at  the  interface  are  perfectly  matched 
at  every  point,  and  if  the  displacement  field  at  the  interface  was  exact,  then  the 
integrand  in  Equation  2.3  is  zero  at  every  instant.  However,  approximations  in  the 
interfacial  displacement  field  due  to  the  parameterization  discussed  earlier,  or  due  to 
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a  mismatched  mesh  (leading  to  mismatched  displacement  fields  on  the  two  surfaces 
forming  the  interface)  will  cause  the  integrand  to  be  non-zero. 

In  the  methodology  developed,  the  decomposed  solution  is  obtained  by 
minimizing  the  accumulated  error  of  Equation  2.3.  Thus  for  a  load  step  over  time  t, 
the  error  to  be  minimized  becomes: 


F  =  W  +W  +W  +W 

^  vrcon,p  T  vv  Mr  I  comp  T  vr  PCB  T  vrstdrlPCB 

“  XUC  -  f'hr  Uc  +  X  X  Fc,  Auc, 


2.4 


0  /=! 
t  n 


1  1  H 
+— ii^K;>un -t,h  u„ +y  y  f„a«„ 

2  P  P  P  ,hp  P  Pi  Pi 


0  /= 1 


Subscript  c  refers  to  the  component  and  subscript  p  refers  to  the  PCB  (Deshpande 
and  Subbarayan,  1999).  By  choosing  an  incrementally  small  time  step  t,  the  time- 
dependent  response  can  also  be  accurately  captured  in  an  incremental  manner, 
analogous  to  the  iterative  schemes  of  nonlinear  finite  element  analysis. 


The  forces  acting  on  any  node  in  the  package  are  the  sum  of  thermal  forces 
and  mechanical  forces: 

f=f„+f„  2.5 

where  the  subscripts  th  and  R  refer  to  the  thermal  forces  and  mechanical  forces 
respectively.  The  thermal  forces  are  the  forces  developed  at  the  retained  nodes  of 
superelements  due  to  thermal  strain.  The  forces  fR  are  the  result  of  externally  applied 
tractions  at  the  retained  nodes.  The  ABAQUS  finite  element  code  outputs  the  force 
vector  fth  during  the  process  of  building  the  superelements.  Since  thermal  strains  by 
definition  do  not  cause  a  mechanical  stress,  the  contribution  of  the  thermal  forces  is 
subtracted  from  the  stored  energy  in  Equation  2.4. 


35 


If  solder  joints  are  assumed  to  behave  linearly  (for  a  preliminary  analysis), 
then  the  error  to  be  minimized  becomes  the  same  as  the  total  potential  energy  of  the 
system: 


E  =  Wcomp+WPCB  +  WcPi 

-f„  u  „  +-uL  K  u  -f„  u 


xup 


cp.^s^cp,  lth"cPi 


2.6 


As  stated  earlier,  the  component  and  PCB  behavior  can  be  approximated  as  being 
linear.  In  Equation  2.6,  subscript  s  refers  to  the  solder  and  subscript  cp,  of  the  fifth 
and  sixth  terms  refers  to  the  solder  interface  with  component  and  PCB  and,  these 
terms  replace  the  third  and  sixth  terms  of  Equation  2.4.  Vector  ucp  is  an  augmented 
vector  comprising  of  vectors  uc  and  up.  Thus,  for  assumed  displacement  vectors  uc 
and  up  the  error  defined  by  Equation  2.6  can  be  calculated,  and  the  process  iterated 
until  the  displacement  vectors  that  yield  the  minimum  error  are  obtained. 

The  equilibrium  state  of  the  system  as  a  whole  is  approximated  by 
minimizing  the  error  function  of  Equation  2.4  in  the  present  procedure.  The 
evaluation  of  the  terms  in  Equation  2.6  is  very  inexpensive  and  straightforward  once 
the  superelement  stiffness  matrix  is  calculated.  The  computational  expense  in 
evaluation  of  Equation  2.4  is  due  to  the  need  for  evaluation  of  the  work  defined  in 
Equation  2.2  for  the  nonlinearly  behaving  solder  alloy. 


2.3.4  Numerical  Simulations  to  Approximate  Solder  Joint  Response 

Following  the  Simultaneous  Analysis  and  Design  (SAND)  idea  applied  by 
Deshpande  et  al.  (1998),  a  statistically  designed  set  of  simulations  are  carried  out 
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next  to  capture  the  nonlinear  response  of  the  solder  joints  to  changes  in  design  and 
analysis  related  parameters.  For  this  modeled  response  to  be  valid  for  the  widest 
possible  choice  of  future  and  present  packages,  careful  thought  must  be  given  to  the 
input  parameters  and  their  ranges.  These  issues  are  described  in  detail  by  Deshpande 
et  al.  (1998),  and  only  the  points  of  major  difference  with  the  prior  research  are 
described  in  the  brief  outline  below. 

In  the  present  study,  a  range  of  0.0-2.0  microns  for  the  relative  shear 
displacement  ( u )  and  a  range  of  0.5-2.0  microns  for  the  relative  axial  displacement 
(v)  were  chosen.  The  linear  variation  of  the  axial  displacement  represented  by  the 
angle  0was  assumed  to  have  a  range  of  -1°  to  1°.  These  ranges  were  based  on  Moire 
experiments  on  several  225  I/O  PBGA  packages  that  indicated  a  peak  relative  shear 
displacement  of  approximately  1.8  microns,  for  a  temperature  change  of  -80  °C 
(Zhang  et.al,  1998).  In  addition,  work  due  to  local  mismatch  at  each  interface  was 
captured  by  varying  the  top  and  bottom  expansion,  at  and  coefficients  from  7  to 
21  ppm/°C.  A  fractional  factorial  designed  “experiment”  (Montgomery,  1997)  based 
on  the  five  levels  of  these  variables,  as  shown  in  Table  2.1,  was  used  to  capture  the 


response  of  the  solder  joint. 


Table  2.1  The  fractional  factorial  designed  experiment  used  to  build 
regression  model  for  solder  work. 
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Run  u  (gm)  v([im)  e  (radians)  ql  (ppm/°C)  q2  (ppm/°C)  Work  (N-^m) 


1 

2 

2 

0.017455 

21 

21 

25.50 

2 

2 

2 

0 

21 

14 

8.32 

3 

2 

2 

-0.017455 

21 

7 

34.58 

4 

2 

1.25 

0.017455 

14 

21 

24.01 

5 

2 

1.25 

0 

14 

14 

4.92 

6 

2 

1.25 

-0.017455 

14 

7 

33.20 

7 

2 

0.5 

0.017455 

7 

21 

24.33 

8 

2 

0.5 

0 

7 

14 

5.53 

9 

2 

0.5 

-0.017455 

7 

7 

33.65 

10 

1 

2 

0.017455 

14 

14 

26.88 

11 

1 

2 

0 

14 

14 

6.31 

12 

1 

2 

-0.017455 

14 

14 

31.38 

13 

1 

1.25 

0.017455 

14 

14 

25.50 

14 

1 

1.25 

0 

14 

14 

2.30 

15 

1 

1.25 

-0.017455 

14 

14 

30.00 

16 

1 

0.5 

0.017455 

14 

14 

25.73 

17 

1 

0.5 

0 

14 

14 

2.71 

18 

1 

0.5 

-0.017455 

14 

14 

30.25 

19 

0 

2 

0.017455 

21 

21 

28.93 

20 

0 

2 

0 

21 

14 

5.69 

21 

0 

2 

-0.017455 

21 

7 

29.13 

22 

0 

1.25 

0.017455 

14 

21 

27.40 

23 

0 

1.25 

0 

14 

14 

1.11 

24 

0 

1.25 

-0.017455 

14 

7 

27.75 

25 

0 

0.5 

0.017455 

7 

21 

27.74 

26 

0 

0.5 

0 

7 

14 

2.07 

27 

0 

0.5 

-0.017455 

7 

7 

28.22 

A  non-linear  elastic-plastic  analysis  was  performed  on  each  of  the  twenty- 
seven  settings  in  the  designed  experiment  using  the  integrated  solder  shape/life 
prediction  method  developed  by  Subbarayan  (1996).  A  typical  finite  element  model 
of  the  solder  joint  generated  by  using  this  methodology  is  shown  in  Figure  2.4.  For 
given  values  of  the  design  variables,  the  program  developed  by  Subbarayan  predicts 
the  shape  of  the  solder  joints  and  the  shape  is  then  automatically  converted  into  a  3- 
D  finite  element  mesh,  which  is  analyzed  using  the  ABAQUS  finite  element  code. 


Figure  2.4  A  typical  finite  element  model  of  the  solder  joint. 


Linear  regression  models  were  built  to  relate  the  work  quantities  defined  in 
Equation  2.4  with  the  five  input  parameters.  These  regression  models  represent  a 
stand-alone  module  for  solder  joints.  Even  though  the  five  inputs  to  this  module  are 
analysis-related  parameters  ( u  and  v  displacements,  0,  ai  and  (X2),  it  is  quite  simple 
to  add  design-related  variables  such  as  the  solder  volume  and/or  solder  pad 
diameters  as  additional  variables  to  this  module.  Such  a  module  can  be  used  in  the 
design  and  analysis  of  packages  with  widely  different  geometries  as  discussed  by 
Deshpande  et.al.  (1998).  In  general,  the  regression  models  are  possible  only  under 
monotonic  loading  conditions  since  the  stress-strain  relation  and  therefore  the 
inelastic  dissipation  is  non-unique  (or  history  dependent)  under  reversed  loading 
conditions.  Thus,  while  the  developed  methodology  is  general  and  applicable  to 
reversed  loading  conditions,  the  regression  model’s  validity  is  limited  to  monotonic 
loading  conditions. 


2.3.5  Code  and  Code  Input  Deck 
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The  code  developed  must  be  general  enough  to  solve  many  types  of 
electronic  packages.  A  major  aspect  of  the  implementation  scheme  is  the  underlying 
optimization  procedure  that  iteratively  determines  the  solution  to  the  unknowns  for 
any  area  array  package.  Our  code  incorporates  the  non-linear  optimization  code 
NLPQL  (Schittkowski,  1985).  This  code  iteratively  improved  the  variables  until  the 
optimum  point  was  reached.  The  NLPQL  routine  finds  the  minimum  or  maximum 
of  a  nonlinear  function  of  n  optimization  variables  subject  to  nonlinear  or  linear 
constraints  using  a  Sequential  Quadratic  Programming  (SQP)  method.  The  gradients 
were  calculated  at  each  iteration  by  forward  differences.  However,  as  the  size  of  the 
packages  increase,  so  will  the  number  of  optimization  variables  and  number  of 
function  evaluations.  Use  of  numerical  gradients  subroutines  with  NLPQL  require 
(n+1)  function  evaluations  at  each  iteration.  For  example,  a  225  I/O  area-array 
package  with  36  solder  joints  at  1/8  scale  (288  optimization  parameters)  would 
require  289  function  evaluations  at  each  iteration.  This  can  be  computational 
expensive  with  CPU  times  approaching  or  exceeding  that  of  a  full  three-dimensional 
non-linear  finite  element  analysis.  In  implementing  our  scheme,  we  developed  an 
analytical  gradient  subroutine  for  the  package’s  linear  substructures.  This  reduced 
the  number  of  function  and  gradient  evaluations  to  one  each  for  a  total  of  two 
evaluations  per  iteration,  greatly  decreasing  CPU  times.  The  gradient  of  the  error 
function  (Equation  2.4)  is  calculated  in  the  code  as: 


dE_ 

dx 


If  .. 


du„ 

K  u  — p- 
p  p  dx 


»  dx 


dWdiy 
du  dx 
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where  W  is  the  regression  model  for  the  accumulated  work  done  by  the  solder  joint. 
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A  typical  input  deck  needed  to  successfully  run  the  code  is  shown  in  Figure 
2.5.  With  this  input  deck,  the  code  generates  the  displacement  field  at  each  interface 
by  first  matching  the  retained  interfacial  nodes  with  the  appropriate  solder  joint 
center  node.  The  program  then  proceeds  to  calculate  the  distance  from  the  package 
neutral  point  to  each  joint  center,  Rnp,  and  the  distance  from  the  matching  nodal 
points  to  the  joint  center,  Rnode-  By  determining  where  the  joint  centers  lie  at 
solder/component  interface  or  solder/PCB  interface,  the  code  generates  the 
displacement  fields.  Figure  2.6  illustrates  the  flow  of  control  within  the  developed 
procedure. 

Typical  Code  Input  Deck 

Location  of  package  neutral  point 
#  of  solder  joints  in  package 
*  Fraction  of  solder  joint  in  package 

Node  #  and  x,y,z  coordinates  of  solder  joint  centers 
Node  #  and  x,y,z  coordinates  of  retained  solder  nodes 
9  Upper  and  lower  solder  joint  pad  diameters 
Height  of  PCB,  component,  and  solder  joint 
9  #  of  retained  nodes  for  PCB  and  component 

9  Node  #  and  x,y,z  coordinates  of  retained  PCB  and  component  nodes 
9  Package  AT 

9  x,  y,  z  thermal  expansion  coefficients  of  component,  solder,  and  PCB 
Super  Element  stiffness  matrices  for  PCB  and  component 
Equations  to  calculate  solder  work  for  nonlinear  solder 
Thermal  force  data  for  PCB  and  component 


Figure  2.5  Typical  code  input  deck. 
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Figure  2.6  Flow  control  in  the  program. 

2.3.6  Solving  for  the  Displacements 

The  decomposition  methodology  described  earlier  was  applied  systematically 
to  two  2-D  electronic  packages,  Deshpande  and  Subbarayan’s  3-D  5x5  hypothetical 
electronic  package,  and  a  225  I/O  PBGA  electronic  package  as  described  below.  For 
the  purposes  of  verification  of  the  code,  the  solder  material  was  initially  assumed  to 
behave  linearly.  This  is  not  a  limitation  of  the  code,  but  an  assumption  that  enabled 
quick  testing  of  the  developed  code.  The  finite  element  models  of  the  2-D  packages 
are  shown  in  Figure  2.7.  The  first  package  has  only  one  solder  joint  (2  interfaces 
with  8  unknown  variables)  and  the  second  package  has  3  solder  joints  (6  interfaces 
with  24  unknown  variables).  The  full  three-dimensional  finite  element  model  of  the 
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5x5  hypothetical  area-array  package  is  shown  in  Figure  2.8.  The  cross-section  of 
this  package  is  shown  in  Figure  2.9.  Because  of  symmetry,  only  one-eighth  of  the 
package  was  modeled  which  had  6  solder  joints  (12  interfaces  with  48  unknown 
variables).  A  full  three-dimensional  finite  element  model  of  the  225  I/O  PBGA 
package  was  also  built  and  is  shown  in  Figure  2.10.  The  cross-section  of  this 
package  is  shown  in  Figure  2.11.  A  one-eighth  symmetry  model  of  the  225  I/O 
package  has  36  solder  joints  (72  interfaces  with  288  unknown  variables). 


2-D  Single-Joint 


2-D  Three -Joints 


Figure  2.7  Finite  element  models  of  hypothetical  2-D  electronic 
packages. 


Figure  2.8  The  3-D  FE  model  of  the  5x5  package. 


Figure  2.9  The  cross-section  of  the  hypothetical  5x5  model. 
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Figure  2.10  The  3-D  FE  model  of  the  225  I/O  PBGA  electronic  package. 


Figure  2.11  Cross-section  of  the  225  I/O  PBGA  electronic  package. 

The  material  properties  for  the  2-D  models  are  shown  in  Table  2.2. 
The  material  properties  for  the  3-D  models  were  chosen  to  be  the  same  as  that  given 
in  Deshpande  et.al.  (1998)  and  are  shown  in  Table  2.3.  The  2-D  models  were 
subjected  to  a  uniform  temperature  rise  of  100°C  from  the  initial  temperature  of  0°C 
and  the  3-D  models  were  subjected  to  an  80°C  temperature  rise  from  the  same  initial 


temperature. 


Table  2.2  Material  properties  used  in  2-D  electronic  packages. 


Material 

EfSISBHEBffli 

CTE  (In-Plane) 

PCB  (FR-4) 

20 

15.5 

0.20 

1.47E-05 

5.12E-05 

Component  (Silicon) 

20 

131 

0.28 

2.60E-06 

2.60E-06 

Solder  <63Sn37Pb) 

0 

32 

0.35 

2.10E-05 

2.10E-05 

20 

30 

50 

27 

100 

22 
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Table  2.3  Material  properties  used  in  3-D  electronic  packages. 


Solder  (63Sn37Pb) 

E  (GPA) 

Poisson’s  Ration 

Temperature  °K 

Yield  Stress  (GPA) 

Plastic  Strain 

Temperature  °K 

34.055 

0.3 

273 

0.0138 

0 

273 

33.329 

0.3 

278 

0.0276 

7.30E-04 

273 

32.603 

0.3 

283 

0.0414 

7.02E-03 

273 

31.877 

0.3 

288 

0.0552 

3.49E-02 

273 

31.151 

0.3 

293 

0.0689 

1.21  E-01 

273 

30.425 

0.3 

298 

0.0827 

3.36E-01 

273 

29.699 

0.3 

303 

0.0965 

7.93E-01 

273 

28.973 

0.3 

308 

0.1100 

1.67E+00 

273 

28.247 

0.3 

313 

0.1240 

3.22E+00 

273 

27.522 

0.3 

318 

0.1380 

5.80E+00 

273 

26.796 

0.3 

323 

0.0138 

0 

373 

26.07 

0.3 

328 

0.0276 

1.62E-02 

373 

25.344 

0.3 

333 

0.0414 

1.56E-01 

373 

24.618 

0.3 

338 

0.0552 

7.76E-01 

373 

23.892 

0.3 

343 

0.0689 

2.70E+00 

373 

23.166 

0.3 

348 

0.0827 

7.46E+00 

373 

22.44 

0.3 

353 

0.0965 

1.76E+01 

373 

21.714 

0.3 

358 

0.1100 

3.71E+01 

373 

20.988 

0.3 

363 

0.1240 

7.16E+01 

373 

20.262 

0.3 

368 

0.1380 

1.29E+02 

373 

19.536 

0.3 

373 

CTE  =  21.0E-6 

PCB 

Die 

Over  Mold 

E  (GPA) 

E  (GPA) 

E  (GPA) 

E  (GPA) 

19.0 

18.2 

130.0 

15.5 

Poisson’s  Ratio 

Poisson’s  Ratio 

Poisson’s  Ratio 

Poisson’s  Ratio 

0.2 

0.19 

0.28 

0.25 

CTE 

CTE 

CTE 

CTE 

1.50E-05 

1.60E-05 

2.62E-06 

1.50E-05 

The  displacements  obtained  with  decomposed  analysis  procedure  were 
compared  with  those  obtained  using  non-decomposed,  full  finite  element  models. 
The  analysis  was  carried  out  using  the  ABAQUS  finite  element  code  (HKS,  1997). 
In  addition,  the  assumption  of  linear  elastic  behavior  in  our  preliminary  analyses, 
gave  us  the  opportunity  to  test  the  non-linear  optimization  code  NLPQL  by  assuming 
all  nodal  displacements  as  unknown  variables  and  comparing  the  results  with  the 
full  finite  element  analysis.  In  this  case,  the  present  procedure  reduces  to  that  of  the 
classical  substructuring  approach  of  finite  elements,  only  the  solution  is  obtained 
using  an  iterative  approach  as  opposed  to  using  matrix  factorization  methods.  The 
error  between  the  iterative  and  non-iterative  solutions  in  these  cases  was  negligible. 
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2.3.7  Error  Measures 


Finally,  a  discussion  of  the  error  measures  used  to  evaluate  the  accuracy  of 
the  decomposed  analysis  procedure  is  in  order.  In  keeping  with  the  parameterization 
scheme  discussed  earlier  and  illustrated  in  Figure  2.2,  the  errors  in  displacements 
may  be  divided  into  two  kinds.  The  first  kind  is  the  error  in  displacement  at  any 
interfacial  node,  expressed  in  a  suitable  non-dimensional  form.  This  is  referred  to  as 
the  displacement  error  or  simply  the  error  in  the  present  paper.  Mathematically,  the 
chosen  error  is  Euclidean  in  form  and  is  defined  as: 


j  M 

AverageError  =  —  \  Error M  2.9 

where  M  is  the  total  number  of  displacement  DOFs  for  component  and  PCB,  N  is  the 
total  number  of  solder  nodes,  y,  is  the  displacement  predicted  by  the  full  FE 
simulation  (the  reference  for  evaluating  the  accuracy  of  the  decomposed  procedure) 
at  the  interface,  and  y,  is  the  solution  from  the  decomposed  analysis  at  the  exact  same 
location. 

The  second  kind  of  error  is  that  in  the  difference  between  displacements  at 
two  points.  This  is  useful  in  calculating  the  relative  shear  and  axial  displacements  of 
the  solder  joints,  which  in  turn  are  important  for  predicting  the  fatigue  life  of  the 
joint  under  cyclic  loading.  This  error  is  denoted  as  the  relative  error  and  its  average 
value  is  mathematically  defined  as: 
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Average  Relative  Error  = 


(  i  \l/2 

2>  relative  V  relative  )  ' 
i= 1 


(  j  v/2 

^  (y relative  ) 


2.10 


V'='  ; 

where  7  is  the  total  number  of  solder  joints,  y' relative  is  the  difference  in  displacements 
between  two  nodes  in  the  full  FE  model,  and  yl  relative  is  the  same  difference 
calculated  using  the  decomposed  analysis.  Since  the  relative  error  measures  the 
error  in  displacement  difference,  it  represents  a  higher  order  error  compared  to  the 
displacement  error. 


2.4  Results  and  Discussion 

Tables  2.4  lists  the  results  of  the  analysis  using  the  decomposed  solution  with 
displacement  fields  as  compared  to  the  finite  element  solution.  The  average 
displacement  error  is  0.74%  for  the  single-joint  2-D  model  and  5.0%  for  the  three- 
joint  2-D  model.  The  average  displacement  error  for  the  3-D  5x5  area-array  model 
with  linear  solder  joints  was  6.04%  while  the  average  relative  shear  and  axial 
displacement  errors  were  6.65%  and  6.03%  respectively.  The  3-D  5x5  area-array 
model  with  nonlinear  solder  joints  had  an  average  displacement  error  of  4.98%  while 
the  relative  shear  and  axial  errors  were  9.00%  and  16.43%  respectively.  This 
illustrates  the  difficulty  in  achieving  small  relative  errors.  Lastly,  the  225  I/O  PBGA 
model  had  an  average  displacement  error  of  5.46%,  while  the  relative  shear  and  axial 
errors  of  41.93%  and  14.42%  respectively.  Decomposed  and  the  full  FE  analysis 
displacement  plots  for  the  5x5  and  225  I/O  arrays  are  shown  in  Figures  2.12  and 


2.13. 
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Table  2.4  Displacement  errors  at  component  and  PWB  interfaces. 


Decomposed  Model 

Displacement  Error  (%) 
u  v  w 

Relative  Displacement 
Error  (%) 

shear  axial 

R2 

U  V  w 

2-D  Single  Joint  -  Component 

2-D  Single  Joint  -  PWB 

2.26  1.34 

0.9997  0.9998 

0.9997  0.9998 

2-D  Three  Joint  -  Component 

2-D  Three  Joint  -  PWB 

9.31  5.01 

0.9855  0.9998 

0.9264  0.9998 

5x5  Array  w/  Linear  Solder  -  Component 

5x5  Array  w /  Linear  Solder  -  PWB 

6.65  6.03 

0.9815  0.9913  0.8135 
0.9837  0.9949  0.8650 

5x5  Array  w /  Nonlinear  Solder  -Component 

5x5  Array  w/  Nonlinear  Solder  -  PWB 

7.63  5.60  5.42 

1.95  2.12  7.12 

Average  4.98 

9.00  16.43 

0.9626  0.9929  0.5263 
0.9974  0.9989  0.3981 

225  PBGA  w/  Nonlinear  Solder-  Component 

225  PBGA  w/  Nonlinear  Solder  -  PWB 

2.29  3.71  8.75 

4.48  3.96  9.58 

Average  5.46 

41.93  14.42 

0.9970  0.9964  0.8694 

0.9814  0.9936  0.8612 

Component  Displacement  Comparison  -  5x5  Array 


♦  u-FE  — u-Decom  posed  -a—  v-FE 

— -a —  ^Decomposed  — * — w-FE  — • — w-Decomposed 


Figure  2.12  Component  displacement  comparison  for  5x5  area-array 
model  depicting  solder  joints  1-6  (J1  -  J6). 
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Component  Displacement  Comparison  -  225  PBGA 


1  3  5  7  9  11  13  15  17  19  21  23  25  27  29  31  33  35  37  39  41  43  45 

Diagonal  Nodes 


u-FE 

— u-Decomposed 

-B-v-FE 

— #—  v- Decomposed 

— ■— w-FE 

— w-Decomposed 

Figure  2.13  Component  displacement  comparison  for  225  I/O  PBGA 
Package  depicting  diagonal  solder  joints  1-8  (J1  -  J8). 

For  the  2-D  models,  the  major  effort  in  creating  a  new  code  input  deck  was 
the  generation  of  new  superelement  matrices  for  the  component  and  PCB.  This 
merely  involved  modifying  the  ABAQUS  input  file  to  retain  the  additional  nodes, 
which  took  a  few  minutes.  Each  solder  joint,  although  an  individual  subsystem,  is  a 
“repeated”  structure  within  the  electronic  package,  thus  the  same  superelement 
matrix  can  be  used  for  each  solder  joint.  Exploiting  the  nature  of  these  repeated 
structures  enables  a  considerable  computational  advantage.  Therefore  we  see  by 
simply  changing  the  code  input  deck,  we  are  able  to  carry  out  a  second  analysis. 

As  for  the  3-D  5x5  package,  superelements  were  already  available.  The  code 
was  easily  modified  to  handle  three-dimensional  packages,  as  was  the  code  input 
deck.  For  the  case  with  the  nonlinearly  behaving  solder,  the  regression  model  for  the 
work  defined  in  Equation  2.4  was: 
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W  =  3.05369  +  1.73768  u  +  1.73116  v  -  2.25479  6 

-0.09433  a.  -0.13001a,  -0.13157  u  vA2  -  1.14245  u  92 

1  2  2.11 

- 1. 12784  v  6 2  +0.55165  u 2  -  2.21 197  ud  + 1.36305  vA2 

+  24.235  62 

The  regression  model  for  the  experiment  was  accomplished  using  the  commercially 
available  software  DOE  KISS  (Design-of-Experiment  Keep-It-Simple-Statistically) 
for  Microsoft  Excel  (Digital  Computations,  1997).  The  factors  in  the  equation  are 
represented  by  coded  variables  where  the  low  and  high  levels  of  each  factor  are 
assigned  the  values  of -1  and  +1. 

Table  2.5  shows  the  CPU  times  for  the  decomposed  and  FE  analyses.  The 
CPU  time  required  for  the  3-D  ABAQUS  non-linear  analysis  of  the  5x5  hypothetical 
package  was  5846  seconds  on  a  Pentium-II  400  MHz  dual  processor  machine, 
whereas  the  decomposed  procedure  required  only  37  seconds  of  CPU  time  on  the 
same  machine.  Including  the  CPU  time  for  the  27  experiments  (3923  seconds)  and 
superelements  (4  seconds),  a  47%  savings  in  computational  effort  was  achieved 
through  decomposition.  However,  we  note  that  for  additional  package  constructions 
using  the  same  solder  joint,  the  CPU  time  for  the  27  experiments  is  redundant  since 
the  regression  model  is  reusable!  This  illustrates  the  power  of  the  decomposed 
analysis  approach. 

The  CPU  time  required  for  the  3-D  ABAQUS  non-linear  analysis  of  the  225 
I/O  PBGA  package  was  23782  seconds  on  a  Pentium-13  400  MHz  dual  processor 
machine,  whereas  the  decomposed  procedure  required  only  907  seconds  of  CPU 
time  on  the  same  machine.  Including  the  CPU  time  for  building  the  superelements 
(2195  seconds)  but  ignoring  the  time  required  to  build  the  regression  model  (for 
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reasons  mentioned  above),  a  667%  speedup  in  computational  effort  was  achieved 
through  decomposition. 

Table  2.5  CPU  times  (seconds). 


5x5 

Decomposed 

Model 

225 

Decomposed 

Model 

Super  Elements 

4 

2195 

Solder  DOE 

3923 

0 

Decomposed  Analysis  w/  FD  Gradients 

312 

N/A 

Decomposed  Analysis  w/  Analytic  Gradients 

37 

907 

Total  w/  FD  Gradients 

4239 

N/A 

Total  w /  Analytical  Gradients 

3964 

3102 

3-D  Finite  Element 

5846 

23782 

CPU  Savings 

47% 

667% 

There  are  two  main  reasons  for  these  computational  savings.  The  first  reason 
is  that  the  3-D  models  have  a  much  larger  number  of  nodes  and  elements 
necessitated  by  the  need  to  match  the  refined  solder  mesh  in  the  PCB  and  the 
component.  The  full  FE  model  of  the  225  I/O  PBGA  package  had  over  94000  nodes 
and  over  90000  elements,  whereas  in  the  decomposed  analysis,  the  component  and 
PCB  superelements  were  coarsely  meshed  and  together  they  contained  less  than 
11,000  elements  and  12,000  nodes.  The  second  reason  is  that  the  problem  with 
many  unknown  degrees  of  freedom  at  the  interface  has  been  reduced  to  only  four 
unknowns  by  approximating  the  displacement  fields  at  the  interface  through  the 
chosen  parameterization  scheme.  Thus,  the  parameterization  represents  a  set  of  basis 
functions  for  the  displacement  at  the  interface,  chosen  with  the  benefit  of  physical 
intuition. 

The  main  sources  of  errors  in  the  displacements  can  come  from:  statistical 
errors  in  the  regression  model,  finite  element  solution  errors  in  capturing  the  non¬ 
linear  solder  behavior,  and  the  errors  in  the  ability  of  the  chosen  basis  functions  to 
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represent  the  exact  displacement  field  at  the  interface.  The  regression  model  used 
here  had  a  goodness  of  fit  measured  by  an  R  value  of  higher  than  0.99.  This 
suggests  that  the  regression  model  is  unlikely  to  be  a  major  source  of  error.  To 
check  if  material  behavior  could  be  a  source  of  error,  analysis  of  the  5x5  area  array 
package  was  repeated  by  assuming  linear  elastic  material  behavior  for  solder  joints. 
In  this  case,  since  the  behavior  was  linear,  a  superelement  was  built  for  the  solder 
joint  and  the  response  of  the  solder  joint  was  calculated  as  described  by  Equation 
2.6.  Thus  closed  form  equations  for  work  done  were  available  for  all  the  three 
subsystems  and  no  regression  model  was  necessary.  As  stated  earlier,  the  average 
relative  shear  and  axial  errors  were  6.65%  and  6.03%  respectively  (See  Table  2.4). 
However,  the  average  overall  displacement  error  with  linear  solder  behavior  is 
approximately  the  same  as  the  average  error  in  displacements  with  the  non-linear 
solder  behavior,  around  5.0  -  6.0%.  This  suggests  that  an  inaccuracy  in  capturing 
the  non-linear  solder  material  behavior  is  unlikely  to  be  the  major  source  of  error  in 
the  displacements  and  that  the  major  source  of  error  is  due  to  the  approximation 
introduced  by  representing  the  interfacial  displacements  using  linearly  varying  shear 
and  axial  displacements. 

The  225  I/O  PBGA  package  had  comparatively  large  axial  displacement  and 
relative  shear  displacement  errors.  Figure  2.14  (reproduced  from  Zhang  et  al.,  1999) 
shows  the  relative  shear  displacement  along  the  diagonal  joints  of  the  225  I/O  PBGA 
package  determined  experimentally  through  moire  interferometry,  and  compared  to 
axisymmetric,  full  3-D  finite  element,  and  decomposed  models.  The  large  error  in 
the  decomposed  solution  beyond  the  die  region  is  due  to  this  region’s  relative 
contribution  to  the  overall  error  function  of  Equation  2.4.  In  Figure  2.15,  the 
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relative  shear  displacement  errors  is  overlaid  on  the  error  gradient  with  respect  to 
variable  yo  which  accounts  for  the  relative  shear  displacement  of  the  solder  joint  due 
to  thermal  mismatch.  As  seen  in  the  figure,  the  gradients  within  the  die  area  of  the 
package  are  larger  than  the  gradients  outside  of  the  die  area.  This  can  be  attributed 
to  the  die  area  being  stiffer  (die  has  larger  Young’s  modulus)  than  the  rest  of  the 
package.  Since  the  iterative  steps  of  the  NLPQL  routine  depend  on  the  magnitude  of 
the  gradients,  it  is  natural  that  the  larger  gradients  in  the  die  region  receive  a  greater 
attention.  In  Figure  2.15,  it  can  be  seen  that  the  relative  shear  displacement  errors 
for  the  solder  joints  within  the  die  region  of  the  package  are  between  10%  to  15% 
while  the  errors  outside  of  the  die  area  can  range  from  25%  to  85%. 

One  may  attempt  to  scale  the  gradients  in  such  a  manner  that  they  are  all  of 
equal  magnitude,  but  such  a  step  is  unnecessary  if  one  is  interested  in  identifying  the 
critical  solder  joint.  This  is  since,  the  solder  joint  with  maximum  relative 
displacement  is  at  the  die  edge  (joint  number  3  from  center)  and  as  Figure  2.14 
shows  the  decomposed  relative  shear  displacement  at  this  joint  agrees  well  with  the 
full  3-D  finite  element  model  and  experimental  data. 


Relative  Displacement  (microns) 
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Diagonal  Solder  Joint  Number 


Figure  2.14  Relative  shear  displacement  comparisons 
(Zhang,  et  al.,  1999). 
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DNP  (mm) 
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Figure  2.15  225  PBGA  shear  gradient/relative  error  plot. 

Finally,  we  found  that  the  use  of  explicit  gradient  subroutines  decreased  the 
CPU  time  for  the  decomposed  solution  procedure  significantly.  Analysis  of  the  5x5 
area-array  model  using  both  finite  difference  and  analytic  gradients  showed  a  275 
second  decrease  in  CPU  time  from  312  to  37  seconds.  For  a  larger  package,  such  as 
the  225  I/O  PBGA,  the  decomposed  analysis  proved  not  to  be  practical  without  the 
analytical  gradients. 

2.5  Conclusions 

The  present  study  has  demonstrated  that  using  the  decomposition  method  and 
a  more  structured  approach  to  generating  displacement  fields,  an  accurate  solution 
can  be  obtained  with  substantial  savings  in  the  computational  effort.  For  example,  a 
667%  speedup  was  achieved  using  the  decomposed  analysis  procedure  on  the  225 
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I/O  PBGA  package  compared  to  a  three-dimensional  nonlinear  finite  element 
analysis.  Deshpande  and  Subbarayan  (1999)  showed  the  methodology  to  be  scalable 
to  arbitrary  package  sizes  with  potentially  orders  of  magnitude  savings  in 
computational  effort.  The  methodology  can  be  applied  to  a  wide  range  of  package 
constructions  by  just  providing  a  new  code  input  deck.  Currently  research  is 
underway  to  develop  methods  to  further  reduce  the  relative  error. 


Chapter  3 


A  Model  For  Assessing  The  Shape  Of  Solder  Joints  In  The  Presence 
Of  Board  Warpage  And  Volume  Variation  In  Area- Array  Packages 

3.1  Chapter  Overview 

In  this  chapter  we  introduce  novel  methodology  studying  the  effect  of 
printed  circuit  board  (PCB)  warpage  and  solder  volume  variation  on  the  final 
equilibrium  configuration  on  area-array  packages.  The  developed  procedure 
is  demonstrated  on  a  three-dimensional  hypothetical  electronic  package. 

The  objectives  of  this  chapter: 

•  Apply  homogeneous  transformation  theory  to  predict  equilibrium 
configuration  of  electronic  packages. 

•  Use  rigid  rotations  of  package. 

•  Predict  equilibrium  in  presence  of  arbitrary  warpage. 

•  Validate  methodology. 

3.2  Introduction 

Area-array  packages  in  general  are  increasingly  popular  choices  in  the  drive 
towards  cost  reduction  and  miniaturization.  They  offer  enormous  area  reductions  in 
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comparison  to  quad  flat  packages  (QFPs)  and  increasing  potential  to  do  so  without 
adding  to  system-level  costs.  As  with  any  type  of  electronic  package,  reliability  is 
one  of  the  major  concerns,  with  low  cycle  fatigue  of  solder  joints  being  one  of  the 
most  common  failure  mechanisms.  Numerical  techniques  such  as  finite  element 
analysis  are  commonly  used  to  determine  package  reliability.  These  reliability 
assessment  techniques  normally  consist  of  three  basic  steps  to  calculate  (Subbarayan, 
et  al.,  1999) 

1.  solder  joint  shapes, 

2.  thermal  stress/strain  distribution  in  the  most  susceptible  joint  after 
solidification  and  during  thermal  cycling,  and 

3.  fatigue  lives  of  the  most  susceptible  or  any  other  selected  solder  joints. 

The  ultimate  goal  of  the  present  research  is  to  understand  and  systematically 

study  the  reliability  implications  of  component  and  circuit  board  warpage  and  solder 
volume  variations  since  it  is  believed  that  warpage  could  lead  to  non-intuitive  failures 
in  joints  other  than  the  usual  “comer  joint”  at  the  die  edge.  However,  this  chapter 
deals  with  the  first  step  mentioned  above.  The  end  result  is  a  tool  to  accurately 
characterize  the  shapes  of  solder  joints  affected  by  design  specifications,  e.g.,  solder 
height,  volume,  circular  pad  diameter,  weight  of  the  package  borne  by  solder  joint, 
and  manufacturing  variations,  e.g.  warpage. 

Solder  joint  shape  resulting  from  warped  circuit  boards  have  not  received 
much  attention  in  the  literature.  Chan,  et  al.,  (1997)  and  Tower,  et  al.,  (1999)  appear 
to  be  the  only  studies  of  this  nature.  Chan,  et  al.,  established  a  reliability  model  to 
study  thermal  fatigue  in  plastic  ball  grid  array  assemblies.  The  model  was  used  to 
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simulate  a  configuration  with  a  large  number  of  warpage  affected  solder  joints. 
Tower,  et  al.,  developed  a  model  to  predict  yield  of  flip-chip  solder  assemblies  by 
considering  mean  and  standard  deviation  of  volume  distribution,  the  assembly 
warpage,  the  die  size,  and  the  number  of  input/outputs.  These  studies,  however, 
ignored  the  pad  rotation  that  would  naturally  result  from  warped  circuit  boards. 
Renken  and  Subbarayan  (2000)  have  shown  that  a  5°  pad  rotation  could  nearly  halve 
the  standoff  height  in  area-array  solder  joints  due  to  the  resulting  nonsymmetric 
distribution  of  droplet  volume  (See  Figure  3.1).  Thus,  the  inclusion  of  pad  rotation  in 
the  present  study  is  believed  to  have  important  consequences  for  the  solder  joint 
reliability.  The  developed  methodology  is  general  in  nature  and  is  demonstrated  on  a 
hypothetical  3x3  area-array  package  shown  in  Figure  3.2.  The  size  of  the  package  is 
1.2  cm  x  1.2  cm,  and  the  diameter  of  the  solder  pads  on  the  package  and  PCB  is  0.375 
mm.  The  number  of  I/O’s  is  9  with  a  5.0  mm  pitch. 


Height  before  tilt:  57.84  microns 


Figure  3.  2  Hypothetical  3x3  array 
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3.3  Methodology 
3.3.1  Model  Overview 

The  model  developed  in  the  present  study  was  based  on  the  following 
assumptions. 

1 .  Circular  pads  of  the  same  size  were  used  on  both  the  package  and  PCB. 

2.  The  assembly  warpage  value  and  subsequent  warpage  shape,  e.g.,  arch, 
bowl,  or  saddle,  of  the  package  and  PCB  are  known  (See  Figure  3.3). 

3.  Solder  mean  volume  and  standard  deviation,  based  on  a  normal 
distribution,  can  be  given  as  inputs. 

The  static  equilibrium  of  the  system  is  achieved  when  the  forces  and  moments  on  the 
package  are  balanced.  These  forces  and  moments  result  from  the  restoring  forces 
and  moments  of  the  individual  solder  joints. 


Saddle-Type  Warpage 


Figure  3.3  Types  of  warpage. 
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As  illustrated  in  Figure  3.4,  on  a  single  joint,  the  restoring  shear  and  normal 
forces  together  with  the  restoring  moment  balance  the  applied  load,  Fa,  in  the  z 
direction  and  the  applied  moment,  Ma,  about  the  x-axis.  However,  the  package 
rotation/translation  together  with  the  PCB  warpage  will  orient  the  top  and  bottom 
pads  in  a  configuration  that  is  quite  different  from  the  reference  configuration  shown 
in  Figure  3.4.  Thus,  the  package  transformation  and  the  PCB  warpage  have  to  be 
mathematically  described  and  their  effects  transformed  to  the  reference  configuration 
so  the  restoring  forces  and  moments  can  be  calculated.  In  the  following  sections,  the 
warpage,  the  package  motion,  the  calculation  of  restoring  forces  and  the 
transformation  of  warped  configuration  to  the  reference  one  are  described  in  detail. 

Z 

Fs:  Restoring  Shear 
Fn:  Restoring  Normal 
Mr:  Restoring  Moment 
W:  Package  Weight 
Ma:  Applied  Moment 
H:  Height 
9:  Tilt  Angle 

PCB  I 


Figure  3.4  Energy  balance  for  tilted  solder  joint. 
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3.3.2  Description  of  the  Warpage 

The  warpage  of  the  PCB  or  the  package  can  be  mathematically  described  as 
the  implicitly  defined  surface: 

f(X,Y,Z)  =  0  '  3.1 

where,  the  warpage  magnitude  Z  varies  from  one  solder  joint  location  to  another. 
Now,  we  can  assign  a  local  coordinate  system  to  the  bottom/top  pad  of  each  solder 
joint.  Assuming  initially  that  there  is  no  warpage,  and  therefore  the  local  axes  are 
aligned  with  the  global  coordinate  axes,  the  transformation  due  to  the  PCB  warpage 
can  be  described  as  a  rotation  matrix: 

nx  nx  nx  0 

Yly  fly  fly  0 

nz  nyz  nz  0 

0  0  0  1 

Where,  the  non-zero  terms  in  the  first  three  columns  are  the  direction  cosines  of  the 
transformed  local  coordinate  axes  with  respect  to  the  global  axes  denoted  by  the 
subscripts.  Since  the  gradient  of  the  warped  surfaced  defined  by  Equation  3.1  points 
in  a  direction  normal  to  the  surface,  the  following  is  true: 


Now,  the  only  requirement  for  the  two  vectors  nx  and  ny  (defined  similar  to 


nz  above)  is  that  they  lie  in  the  tangent  plane  defined  by  the  set  of  points  s  such  that: 

V/  •  s  =  0  3.4 
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Assuming  without  any  loss  of  generality  that  the  intersection  of  the  above  tangent 
plane  with  a  plane  parallel  to  the  global  X-Z  plane  as  the  direction  of  the  new  local 
x-axis,  we  can  define  nx  as: 

3.5 

Clearly,  with  the  above  definition,  nx  is  orthogonal  to  nz.  Finally,  ny  can  be  obtained 
as  the  cross  product: 

ny=nzxnx  3.6 

If  the  function  had  been  defined  parametrically,  as  is  customary  in  the  field  of 
geometric  modeling  (Mortensen,  1997),  then  the  surface  definition  of  Equation  3.1 
would  be  in  terms  of  two  parameters  u  and  v  in  the  range  [0,1]: 

X(u,v) 

S (u,v)  =  <  Y(u,v)  >  3.7 

Z(u,v ) 

The  vectors  nx,  ny  and  nz  are  easily  obtained  from  this  parametric  definition  as: 

nx=— ,  ny  =  —  and  nz=nxxny  3.8 

du  dv 

3.3.3  Description  of  the  Package  Configuration 

The  package  can  in  general  translate  and  rotate  relative  to  its  original 
configuration  during  the  self-assembly  process.  The  rotation  of  the  package  can  be 
defined  about  an  arbitrary  spatial  axis  originating  at  the  centroid  and  described  by  its 
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three  direction  cosines.  This  rotation  can  be  easily  visualized  as  a  sequence  of 
rotations  about  the  three  coordinate  axes.  The  result  of  the  procedure  is  a  composite 
rotation  matrix: 

8,=8,»^||i|<S;(a)  3.9 

Where  the  R’s  are  the  standard  rotation  matrices  in  homogeneous  coordinates  (Craig, 
1989): 


Rz(<x)  = 


Cos(a )  —  Sin(a )  0  0 
Sin{a)  Cos(a )  0  0 

0  0  10 

0  0  0  1 


*r(P) 


Cos(p)  0  Sin(p)  0 
0  10  0 
-  Sin(p)  0  Cos{p)  0 
0  0  0  1 


10  0  0 

0  Cos(y)  -  Sin(y)  0 
0  Sin(y )  Cos(y)  0 

0  0  0  1 


RA.V)  is  the  transformation  matrix  that  rotates  the  axis  a-deg  about  the  reference  z- 
axis  and  places  the  rotation  axis  in  x-z  plane  of  the  reference  coordinate  system. 
/?y(p)  is  the  transformation  matrix  that  rotates  the  axis  (f-deg  about  the  reference  y- 
axis  and  makes  it  coincident  with  reference  x-axis.  Rx{ y)  is  the  transformation 
matrix  that  rotates  the  package  y-deg  CCW  about  the  reference  x-axis.  The  three 
angles  a,(3,y  are  the  unknowns  to  be  determined  to  describe  the  equilibrium 


orientation  of  the  package. 
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The  translation  of  the  package,  which  in  turn  will  be  imposed  on  the  top  pads 
of  the  solder  joints,  can  be  expressed  as: 


Where,  the  last  column  contains  the  translations  of  the  package  centroid  in  the  three 
coordinate  directions.  These  translations  are  again  unknowns  that  need  to  be 
determined.  Combining  the  translation  and  rotation  of  the  package,  the  composite 
transformation  of  the  package  is  obtained  as: 

T=TnRn  3.10 

c  p  p 

Where,  the  order  of  the  operations  ensures  that  rotation  is  carried  out  prior  to 
translation. 

3.3.4  Restoring  Forces  and  Moments 

Surface  energy  minimization  theory  for  predicting  molten  solder  droplet  shape 
is  well  established  at  the  present  time  (see  Heinrich  (1997)  for  a  survey,  and 
Subbarayan  (1996)  for  a  specific  example).  In  the  present  study,  the  Surface  Evolver 
code  was  used  to  predict  the  shapes  (Brakke,  1994).  Given  a  user-defined  initial 
surface,  the  program  evolves  toward  a  minimum  energy  profile  by  using  a  gradient 
descent  method  on  a  space  of  admissible  surfaces  to  try  to  find  the  local  minimum  of 
the  energy  function.  The  energies  of  relevance  for  a  droplet  are  surface  tension  and 
gravitational  energies.  For  solder  joint  applications  the  gravitational  energy  of  the 
droplet  itself  is  small  compared  to  the  surface  energy  and  can  be  ignored.  Thus,  under 
quasi-static  conditions,  the  restoring  normal  and  shear  forces  generated  by  the  solder 


1  0  0  Tx 
0  10  Ty 
0  0  1  Tz 
0  0  0  1 
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joint  can  be  calculated  as  the  derivatives  of  energy  with  respect  to  a  motion  in  the 
required  directions: 


dE  _  3A 

dP  H  =Const  dP  H  =Const 


3.11 


dE  dA 

-  ~  <T - 

dH  P=Const  dH  P=Const 


3.12 


E  = 


3.13 


where  Fs,  FN  are  the  shear  and  normal  restoring  forces  respectively  (see  Fig.  3.4).  P 
is  the  misalignment,  H  is  the  solder  joint  height,  E  is  the  local  minimum  energy  of 
the  solder  joint,  A  is  the  local  minimum  area,  and  o  is  the  surface  tension  coefficient 
which  is  assumed  to  be  constant  at  every  point  on  the  surface.  In  addition  to  the 
restoring  forces,  a  restoring  moment,  MR,  also  exists: 


,  dE  dA 

Mr  =  — 

*  de  de 


3.14 


For  a  given  set  of  solder  parameters:  pad  dimensions,  solder  volume,  solder 
height,  and  pad  rotation.  Surface  Evolver  determined  the  stable  shape  after  iterative 
computations  and  calculated  the  restoring  forces  and  moments.  The  surface  tension 
coefficient  value  was  350  dyne/cm,  a  value  most  commonly  used  for  eutectic  Sn/Pb 
solder  and  most  other  solder  materials  (Tower,  et  al.,  1999). 


Since  the  package  may  contain  numerous  solder  joints  each  with  its  own  set 
of  parameter  values,  regression  models  were  developed  to  calculate  the  restoring 
forces  and  moments.  The  models  were  developed  using  a  commercially  available 
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multiple  regression-based  analysis  tool  DOE  KISS  (Design-Of-Experiment  Keep-It- 
Simple-Statistically)  for  Microsoft  Excel  (Digital  Computations,  1997).  A  3-level/3 
factor  (33)  full  factorial  designed  experiment  was  used  to  generate  the  force  and 
moment  models.  The  misalignment  was  ignored  in  building  the  regression  model  to 
reduce  the  computational  expense.  While  this  would  lead  to  an  inaccurate  solution,  it 
was  felt  that  this  inaccuracy  would  be  small  for  small  amounts  of  warpage  or  volume 
variations.  The  factors  and  their  levels  are  shown  in  Table  3.1 


Table  3.1  DOE  factors. 


Factor 

Level  2 

IIS5EI 

■I 

mm 

A.  Solder  Volume  (mm)3 

5.15  E-2 

5.90  E-2 

B.  Solder  Height  (mm) 

0.275 

0.400 

C.  Solder  Pad  Tilt,  0  (deg) 

-5 

0 

+5 

The  results  of  the  27  runs  are  shown  in  Table  3.2.  The  force  and  moment  data  were 
generated  during  each  Evolver  run.  The  force  and  moment  models  generated  by 
DOE  KISS  give  explicit  relations  (polynomial  form)  between  the  force  and  moment 
values  and  the  design  parameters  of  the  solder  joint.  These  models  are  listed  below: 

F,  =  -  0.55590  *  C  -  0. 1 6027  *  A  *  C  +  0.98538  *  A  *  C 
4  3.15 

+  0.00895 A  *  B  * C  +  0.03992 *  A2  *C- 0.48075  *B2*C 

Fn  =  - 15.896  +  7.33854  *  A-  34.401  *  B  -  4.78442  *  A*  B  - 
0.58688 *  A2  +  21.347 *B2  - 1.01 196 *  C2  + 1.56582 *  A2  *  S  -  3.16 

1.1 8379 *  A*  B2- 0.97894 *  A  * C2  + 1 .58 1 17  *  B  * C2 

M R  =  -  0.02326  *  C  +  0.02584  *  A  *  C  -  0.045 1 3  *  5  *  C 
-0.04717A* 5 *C  + 0.01046*  A2  *C +  0.06021* F2  *C 


3.17 
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The  volume,  height,  and  rotation  values  in  the  equations  above  are  the  non- 
dimensional  or  coded  value  of  the  parameters.  The  R  value  characterizing 
goodness  of  the  fit  is  0.9993  for  the  shear  restoring  force.  It  is  0.9994  for 
normal  restoring  force  and  0.9069  for  the  restoring  moment. 

Table  3.2  DOE  results. 


the 

the 


Experiment 

Factors 

Volume  (cm3)  Heiqht  (cm) 

Tilt  (deq) 

Shear  Force  fdvne) 

Results 

Normal  Force  (dvne)  X-Moment  (dvne-cm) 

1 

0.0000665 

0.0525 

5 

-0.132642 

-27.0194 

-0.006742 

2 

0.0000665 

0.0525 

0 

0.000000 

-27.0667 

0.000000 

3 

0.0000665 

0.0525 

•5 

0.132640 

-27.0194 

0.006742 

4 

0.0000665 

0.04 

5 

-0.776000 

-10.1411 

-0.021588 

5 

0.0000665 

0.04 

0 

0.000000 

-10.1645 

0.000000 

6 

0.0000665 

0.04 

-5 

0.776000 

-10.1411 

0.021588 

7 

0.0000665 

0.0275 

5 

-2.081600 

43.9081 

0.187852 

8 

0.0000665 

0.0275 

0 

0.000000 

50.47 

0.000000 

9 

0.0000665 

0.0275 

-5 

2.081600 

43.9081 

-0.187852 

10 

0.000059 

0.0525 

5 

0.021000 

-28.2353 

-0.001176 

11 

0.000059 

0.0525 

0 

0.000000 

-28.2737 

0.000000 

12 

0.000059 

0.0525 

-5 

0.021000 

-28.2353 

0.001176 

13 

0.000059 

0.04 

5 

0.537000 

-17.1968 

-0.017277 

14 

0.000059 

0.04 

0 

0.000000 

-17.2557 

0.000000 

15 

0.000059 

0.04 

-5 

0.537000 

-17.1936 

0.017277 

16 

0.000059 

0.0275 

5 

-2.071200 

37.746 

0.069080 

17 

0.000059 

0.0275 

0 

0.000000 

39.8438 

0.000000 

18 

0.000059 

0.0275 

-5 

2.071200 

37.746 

-0.069080 

19 

0.0000515 

0.0525 

5 

0.079690 

-28.4622 

0.004286 

20 

0.0000515 

0.0525 

0 

0.000000 

-28.487292 

0.000000 

21 

0.0000515 

0.0525 

-5 

-0.079690 

-28.4622 

-0.004286 

22 

0.0000515 

0.04 

5 

-0.274866 

-23.4961 

-0.010000 

23 

0.0000515 

0.04 

0 

0.000000 

-23.57 

0.000000 

24 

0.0000515 

0.04 

-5 

0.274866 

-23.4961 

0.010000 

25 

0.0000515 

0.0275 

5 

-1 .833460 

25.28346 

0.010200 

26 

0.0000515 

0.0275 

0 

0.000000 

26 

0.000000 

27 

0.0000515 

0.0275 

-5 

1.833460 

25.28346 

-0.010200 

3.3.5  Transformation  to  Reference  Solder  Joint  Configuration 

The  steps  involved  in  transforming  the  solder  droplet  to  the  reference  configuration, 
taking  into  account  the  package  rotation  and  the  PCB  warpage  are  as  follows: 

1.  Transform  the  top  pad  orientation  to  account  for  package  warpage  and 
package  motion: 

R,=TcRwp  3.18 

where  Rwp  is  the  rotation  due  to  warpage  described  in  Equation  3.2  and  Tc  is  as 


described  in  Equation  3.10. 
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2.  Apply  the  following  rotation  to  the  top  and  bottom  pads  so  that  the 
bottom  pad  is  aligned  with  the  global  axes  and  only  the  top  pad  is  rotated. 

R2=RwbRi  3.19 

where  Rwb  is  rotation  due  to  the  warpage  of  the  PCB. 

3.  Rotate  the  droplet  about  the  global  z-axis  until  an  orientation  is  reached 
in  which  the  rotation  of  the  top  pad  may  be  expressed  as  a  rotation  about 
the  global  x-axis.  That  is,  the  same  final  orientation  R3  of  the  top-pad  is 
to  be  achieved  by  the  following  two  equivalent  rotations: 

R3=Rz(0z)R2=Rx(Ox)  3.20 

where 


Cos(07)  - Sin(0z )  0  0 

Sin(0z )  Cos(dz)  0  0 

0  0  10 

0  0  0  1 


10  0  0 
0  Cos(dj  -Sin(6x)  0 
0  Sin(0x )  Cos(0x )  0 

0  0  0  1 


Equating  the  direction  cosine  with  respect  to  the  z-axis  of  each  transformation 
matrix,  that  is,  the  coefficients  of  the  third  columns  of  Rx( 0X)  and  RA 6-JR2,  we 
get  three  consistent  equations  in  the  two  unknowns  6x  and  Qz- 

rI3  Cos(&)  -  r23  Sin(ft)  =  0  3.21 

rn  Sin(#v)  +  r23  Cos(0z)  =  -Sin(6>()  3.22 


r33  =  Cos(0x) 


3.23 


71 


Where,  r,y  are  the  ijth  terms  of  the  R2  matrix.  The  orthonormality  of  the  R2  matrix 
ensures  the  consistency  of  the  above  three  equations.  Thus,  the  rotation  about  the 
x-axis  to  be  input  into  the  regression  model  is  simply, 

6X  =  CosV  r33)  3.24 

3.3.6  Warpage  Types 

There  are  many  choices  for  the  nature  of  warpage  in  PCBs  that  can  be  used 
for  the  demonstration  of  the  developed  methodology.  In  the  present  study, 
experimental  observations  from  literature  were  used  to  guide  warpage  pattern 
selection.  Amagai  observed  warpage  of  a  Micro-Star™  (p*)  BGA  CSP  post¬ 
assembly  at  25  °C  (Amagai,  1999).  He  then  correlated  the  magnitude  of  assembly 
warpage  to  solder  joint  life  through  Finite  Element  Analysis  (FEA).  The  correlation 
was  strong.  Of  significance  to  the  present  study,  the  warpage  of  the  PCB  post¬ 
assembly  was  arch-type  and  was  in  the  9  to  15  pm  range. 

However,  in  the  present  study,  the  focus  is  restricted  to  the  co-planarity  of  the 
package/component  and  PCB  during  the  reflow  process.  The  reflow  temperature  of 
eutectic  solder  is  around  220°C.  Figure  3.5  shows  a  relatively  flat  package  at  this 
temperature,  hence,  no  warpage  is  present.  Experiments  conducted  on  PCBs  by 
Sutherline,  et  al.  (1998),  showed  a  bowl-type  maximum  out-of-plane  displacement 
difference  of  0.629  mm  for  a  203.2  mm  x  177.8  mm  x  1.524  mm  thick  PCB  with 
traces  at  220°C  using  infrared  reflow.  The  combined  non-zero  flatness  of  the 
package  and  PCB  are  shown  in  Figure  3.6  for  which  we  assumed  a  25  pm  bowl¬ 
shaped  warpage  of  the  PCB.  Non-symmetric  warpage  was  modeled  by  allowing  the 
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right  or  left  side  warpage,  as  measured  from  the  package  centerline,  to  vary  from  0  to 
25  pm  (see  Fig.  3.6). 


180  160  140  120  100  80  60  40  20 

Temperature  (C) 

Figure  3,5  Package  warpage  for  “p*BGA”  CSP  (Amagai,  1999) 


Figure  3.6  Combined  warpage  of  package  and  PWB  during 
reflow. 
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3.3.7  The  Optimization  Algorithm 

The  developed  methodology  requires  six  optimization  parameters  to  describe 
package  rotation  about  its  rotation  axis  (3  parameters)  and  package  translation  (3 
parameters).  Here,  since  the  regression  model  was  simplified  by  ignoring 
misalignment,  fewer  parameters  were  required.  Thus,  the  package  equilibrium 
optimization  problem  may  be  formally  posed  as: 

Find  the  equilibrium  height,  h,  tilt  axis,  and  degree  of  tilt  of  the  package / (h,a,P,Y) 
so  that: 

^package  —  +  ^Si )  =  ®  3.25 

i 

X[MR,+AX,x(FN1+Fs,)]  =  0  3.26 

i 

where  F package  is  the  externally  applied  force  on  the  package  (including  its  weight), 
and  F Ni  and  FS(  are  the  normal  and  shear  restoring  force  vectors  at  joint  i.  AX,  is  the 
distance  from  the  solder  joint  to  the  point  about  which  the  moment  is  calculated.  The 
above  two  vector  equations  are  equivalent  to  a  total  of  six  scalar  equations.  The 
nonlinear  equation  solution  was  reformulated  in  the  present  study  as  an  optimization 
problem  so  an  approximate  solution  that  may  not  satisfy  the  equations  exactly  is  also 
allowed.  The  objective  in  the  optimization  is  to  minimize  the  least  squared  error  in 
the  above  equations.  The  solution  was  determined,  starting  from  a  level  package 
configuration  using  the  NLPQL  code  (Schittkowski,  1985).  This  code  iteratively 
improved  the  design  until  the  optimum  point  was  reached.  The  NLPQL  routine 
finds  the  minimum  or  maximum  of  a  nonlinear  function  of  n  optimization 
parameters  subject  to  nonlinear  or  linear  constraints  using  a  Sequential  Quadratic 
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Programming  (SQP)  method.  The  gradients  were  calculated  at  each  iteration  by 
forward  differences.  For  an  assembly  with  a  given  warpage  level,  the  algorithm  to 
calculate  the  heights  of  all  the  solder  joints  and  package  tilt  is  shown  in  the 
flowchart.  The  height  of  the  solder  joints  due  to  warpage  is  calculated  using 

Z(X,  Y)  =  Zeenter  -  AZ  (X2  +  Y2)/Lpckg  3.27 

The  warpage  is  characterized  by  the  vertical  height  variation  (AZ)  across  the 
diagonal  length  (LpCkg)  of  the  package.  Figure  3.7  illustrates  the  flow  of  control 
within  the  developed  procedure. 


Guess  the 

height  of  the  package 


Guess  the  tilt  axis, 
degree  rotation,  and 

I  translations  of  the  package 


Calculate  the  heights  of 
all  solder  joints  due  to 
warpage  and  package  tilt 


Use  regression  models  to 
calculate  restoring  forces 

and  moments 


Calculate  relative  tilt  and 
offset  of  each  solder  joint 

I  aligned  with  DOE  tilt  axis 


Figure  3.7  Flow  of  control  in  program. 


3.3.8  Validation  of  the  Model 

The  hypothetical  3x3  area-array  electronic  package  was  used  to  test  the  code. 
Eight  test  cases  were  selected  based  on  variations  in  volume  and  warpage: 

1.  Constant  volume  -  no  warpage 

2.  Constant  volume  -  25  pm  left  warpage/0  pm  right  warpage 
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3.  Constant  volume  -  0  (im  left  warpage/25  |im  right  warpage 

4.  Constant  volume  -  25  |im  left  warpage/25  |im  right  warpage 

5.  Volume  varies  linearly  from  left  to  right  -  no  warpage  -> 

6.  Volume  varies  linearly  from  top  to  bottom  -  no  warpage  V 

7.  Volume  varies  linearly  from  left  top  comer  to  right  bottom  comer  -  no 
warpage  ^ 

8.  Volume  varies  linearly  from  left  bottom  comer  to  right  top  comer  -  no 
warpage  7! 

The  results  of  the  test  cases  are  shown  in  Table  3.3  and  Figures  3.8-3.15.  The 
equilibrium  configuration  for  each  test  was  exactly  as  predicted. 


Table  3.3  Test  results  (V  =  Mean  Volume  (5.90E-5  cm3)). 


Volume  (cm3) 

Final  Heiqht  (cm) 

Final  Tilt  (deq)  1 

1 .00*V 

1.00*V 

1.00*V 

0.035 

0.035 

0.035 

0 

0 

0 

Test  #1 

1 .00*V 

1  .oo*v 

1  .oo*v 

0.035 

0.035 

0.035 

0 

0 

0 

i.oo*v 

1.00*V 

1  .oo*v 

0.035 

0.035 

0.035 

0 

0 

0 

i  .ocrv 

1  .oo*v 

1  .oo*v 

0.0343 

0.0356 

0.0347 

0.3314 

0.1252 

0.1200 

Test  #2 

i.oo*v 

1  .oo*v 

1  .oo*v 

0.0355 

0.0357 

0.0347 

0.1666 

0.1200 

0.1200 

i.ocrv 

1  .oo*v 

1  .oo*v 

0.0343 

0.0356 

0.0347 

0.3314 

0.1252 

0.1200 

i  .oo*v 

1  .oo*v 

1  .oo*v 

0.0347 

0.0356 

0.0343 

0.1200 

0.1252 

0.3314 

Test  #3 

1  .oo*v 

1  .oo*v 

1  .oo*v 

0.0347 

0.0357 

0.0355 

0.1200 

0.1200 

0.1666 

1  .oo*v 

1.00*V 

1  .oo*v 

0.0347 

0.0356 

0.0343 

0.1200 

0.1252 

0.3314 

1  .oo*v 

1.00*V 

1  .oo*v 

0.0342 

0.0354 

0.0342 

0.4051 

0.2865 

0.4051 

Test  #4 

1  .oo*v 

1.00*V 

1.00*V 

0.0354 

0.0367 

0.0354 

0.2865 

0.0004 

0.2865 

1.00*V 

1  .oo*v 

1  .oo*v 

0.0342 

0.0354 

0.0342 

0.4051 

0.2865 

0.4051 

1 .04*V 

1.00*V 

0.96*V 

0.0356 

0.0350 

0.0343 

0.0795 

0.0795 

0.0795 

Test  #5 

1 .04*V 

1  .oo*v 

0.96*V 

0.0356 

0.0350 

0.0343 

0.0795 

0.0795 

0.0795 

1 .04*V 

ixxrv 

0.96*V 

0.0356 

0.0350 

0.0343 

0.0795 

0.0795 

0.0795 

1 .04*V 

1 .04*V 

1 .04*V 

0.0356 

0.0356 

0.0356 

0.0795 

0.0795 

0.0795 

Test  #6 

1 .00*V 

1.00*V 

1.00*V 

0.0350 

0.0350 

0.0350 

0.0795 

0.0795 

0.0795 

0.96*V 

0.96*V 

0.96*V 

0.0343 

0.0343 

0.0343 

0.0795 

0.0795 

0.0795 

1 .04*V 

1 .02*V 

1.00*V 

0.0357 

0.0353 

0.0350 

0.0562 

0.0562 

0.0562 

Test  #7 

1 .02*V 

1  .oo*v 

0.98*V 

0.0353 

0.0350 

0.0346 

0.0562 

0.0562 

0.0562 

1.00*V 

0.98*V 

0.96*V 

0.0350 

0.0346 

0.0343 

0.0562 

0.0562 

0.0562 

1 .00*V 

0.98*V 

0.96*V 

0.0350 

0.0346 

0.0343 

0.0562 

0.0562 

0.0562 

Test  #8 

1 .02*V 

1  .oo*v 

0.98*V 

0.0353 

0.0350 

0.0346 

0.0562 

0.0562 

0.0562 

1 .04*V 

1 .02*V 

1.00*V 

0.0357 

0.0353 

0.0350 

0.0562 

0.0562 

0.0562 

Tilt  Axis 


Figure  3.8  Test  #1  results  showing 
tilt  axis  and  degree  tilt. 


3x3  Area-Array  Test  #2 
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Figure  3.9  Test  #2  results  showing 
tilt  axis  and  degree  tilt. 


3x3  Area- Array  Test  #3 
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Figure  3.10  Test  #3  results  showing 
tilt  axis  and  degree  tilt. 


3x3  Area-Array  Test  #4 
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Figure  3.11  Test  #4  results  showing 
tilt  axis  and  degree  tilt. 
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Figure  3.12  Test  #5  results  showing 
tilt  axis  and  degree  tilt. 


Figure  3.13  Test  #6  results  showing 
tilt  axis  and  degree  tilt. 
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3.4  Conclusions 

The  research  has  demonstrated  a  methodology  for  predicting  the  final 
equilibrium  configuration  for  area-array  packages  with  warped  solder  joints  during 
solder  reflow.  This  technique  is  novel  in  its  application  of  homogeneous 
transformation  theory  to  electronic  packages  and  the  use  package  rotation  as  an 
optimization  parameter  to  determine  package  equilibrium  in  the  presence  of 
symmetric  and  non-symmetric  warpage. 

The  methodology  has  been  modified  to  include  package  warpage  and  solder 
joint  misalignment.  The  application  of  the  methodology  to  determine  the  thermal 
stress/strain  distribution  in  the  most  susceptible  joint,  and  the  fatigue  life  of  the  most 
susceptible  or  any  other  selected  solder  joints  in  a  “hypothetical”  electronic  package 
will  be  demonstrated  in  Chapter  4. 


Chapter  4 


A  Methodology  For  Assessing  The  Reliability  Of  Warped  Solder 
Joints  In  Area-Array  Packages 

4.1  Chapter  Overview 

In  this  chapter  we  introduce  a  novel  methodology  for  studying  the 
effect  of  warpage  on  solder  joint  reliability.  The  methodology  consists  of 
combining  the  solder  shape  prediction  methodology  introduced  in  Chapter  3 
with  the  decomposed  analysis  technique  discussed  in  Chapter  2  to  determine 
the  equilibrium  configuration  of  an  area-array  package  with  warped  solder 
joints  and  assessing  the  reliability  of  these  warped  joints. 

Finally,  to  illustrate  the  usefulness,  flexibility  and  strength  of  this 
procedure,  the  procedure  is  demonstrated  on  a  hypothetical  three-dimensional 
electronic  package  with  four  solder  joints. 

The  objectives  of  this  chapter: 

•  Predict  shape  of  solder  joints  in  hypothetical  area-array  package 
using  the  shape  prediction  program  discussed  in  Chapter  3. 

•  Apply  decomposition  technique  to  analyze  a  hypothetical  area- 
array  package  with  warped  solder  joints. 
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•  Finally,  to  assess  the  reliability  of  warped  solder  joints. 

4.2  Introduction 

Reliability  is  one  of  the  major  concerns  for  ball  grid  array  (BGA)  assemblies. 
One  common  failure  mechanism  in  these  packages  is  low-cycle  fatigue  of  solder 
joints.  The  fatigue  life  of  solder  joints  depends  on  the  stress/strain  in  the  solder 
joints  under  the  user  environment.  The  stresses  in  turn  depend  on  the  shape  of  the 
solder  joints.  In  package  design  and  optimization,  fatigue  analysis  derived  by 
experiment,  such  as  accelerated  temperature  cycling,  is  expensive  and  time 
consuming.  A  numerical  model  can  reduce  the  product  development  cycles. 

Restating  from  Chapter  3,  the  numerical  model  should  consist  of  three 
critical  steps: 

1.  solder  joint  shapes, 

2.  thermal  stress/strain  distribution  in  the  most  susceptible  joint  after 
solidification  and  during  thermal  cycling,  and 

3.  fatigue  lives  of  the  most  susceptible  or  any  other  selected  solder  joints. 

The  solder  shape  prediction  model  and  the  code  demonstrated  in  Chapter  3 

accomplish  the  first  step.  The  thermal  stress/strain  calculation  identifies  the  most 
susceptible  solder  joint  and  estimates  its  damage  level  after  each  temperature  cycle. 
The  accumulated  damage  level  can  be  used  to  predict  the  fatigue  life  by  using 
empirical  models  correlated  to  the  large  body  of  existing  experimental  data. 

There  have  been  many  papers  reporting  calculation  techniques  for  the  above 
three  steps.  The  research  in  this  chapter  demonstrates  novel  approaches  for  shape, 
stress/strain,  and  fatigue  calculations.  Since  we  have  already  discussed  the 
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methodology  for  solder  shape  prediction  in  Chapter  3,  this  chapter  will  only  consider 
stress/strain  and  fatigue  calculations. 

4.3  Methodology 

4.3.1  Thermal  Stress/Strain  Models 

Thermal  stress/strain  models  have  been  well  developed  for  surface  mount 
technologies.  It  is  common  to  use  2-D  finite  element  modeling  to  simulate  the 
structures  of  surface  mount  technologies  (Barker,  et  al.,  1993).  For  more  complex 
BGA  assemblies,  3-D  finite  element  modeling  is  necessary  to  simulate  every  solder 
joint  in  order  to  identify  the  joint  most  susceptible  to  failure.  However,  full-scale  3- 
D  modeling  in  not  cost  effective  and  is  usually  replaced  by  alternative  modeling 
techniques  like:  1)  analytical  model  (Borgessen,  et  al.,  1992),  2)  3-D  sliced  FE 
model  (Nagarand  and  Mahakingam,  1993),  3)  macro  (global)/micro  (local)  FE 
model  (Corbin,  1994,  Cheng,  et  al,  1998,  and  Saito,  et  al.,  1999),  4)  Nested  FE 
model  (Darbha  and  Dasgupta,  1999),  and  5)  Decomposition  model  (Deshpande  and 
Subbarayan,  2000).  Zhang,  et  al.,  (1999)  present  an  evaluation  of  the  accuracy  of  2- 
D  and  3-D  structural  approximations  (relative  to  experimental  measurement)  in  the 
analysis  of  area-array  packages. 

In  general,  most  existing  analytical  models  cannot  handle  complicated,  time- 
dependent  deformation  mechanisms  such  as  creep.  The  3-D  sliced  model  cannot 
simulate  all  the  solder  joints.  One  common  implementation  of  the  macro 
(global )/micro  (local)  strategy  replaces  solder  joints  by  equivalent,  nonlinearly 
behaving  beams  in  which  the  beam  properties  are  calculated  using  a  designed 
experiment  based  on  the  analysis  of  a  three-dimensional  solder  model.  Once  the 
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equivalent  beam  has  been  formulated,  a  full-scale  3-D  FE  analysis  of  the  entire 
assembly  is  performed  identifying  the  most  susceptible  solder  joint.  A  micro/local 
model  then  simulates  the  structure  of  that  solder  joint  using  a  much  finer  mesh  than 
the  macro  (global)  model.  The  disadvantage  to  this  approach  is  that  any  component 
or  PCB  change  would  require  a  complete  re-analysis  of  the  3-D  model  with  the 
equivalent  beams.  The  nested  FE  model  (NFEM)  uses  colonies  of  nested  sub¬ 
elements  that  are  created  inside  a  main  element  to  capture  localized  large  gradients 
of  stress  and  strain.  Once  again,  any  change  to  the  assembly  would  require  re¬ 
analysis  of  a  full  3-D  FE  mesh.  The  decomposition  model,  discussed  in  Chapter  2,  is 
free  from  the  macro  (global)/micro  (local)  and  NFEM  disadvantages  since  the 
component  and  PCB  are  separate  subsystems  and  any  changes  would  simply  require 
the  rebuilding  of  that  subsystem  model  alone.  However,  none  of  the  models  above 
appear  to  have  been  used  to  analyze  a  package  with  warped  solder  joints.  The 
decomposition  technique,  with  its  natural  ability  to  accommodate  variations  in 
subsystems  may  allow  one  to  study  the  effects  of  warped  solder  joints  efficiently.  In 
implementing  a  decomposed  analysis  strategy  to  study  warpage  effects  on  solder 
joint  reliability,  the  only  real  need  is  for  a  new  regression  model  for  solder  work  that 
incorporates  design,  analysis  and  manufacturing  parameters.  This  is  since  an 
assumption  is  made  in  the  present  study  that  the  stiffness  matrices  of  PCB  and 
component  are  not  altered  by  their  warpage.  That  is,  the  geometrical  changes  caused 
by  warpage  are  small  enough  that  they  don’t  alter  structural  behavior. 
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4.3.2  Numerical  Simulations  to  Approximate  Solder  Joint  Response 

A  statistically  designed  set  of  simulations  is  carried  out  next  to  capture  the 
nonlinear  response  of  the  solder  joints  to  changes  in  design,  analysis,  and 
manufacturing  related  variables.  For  this  modeled  response  to  be  valid  for  the 
widest  possible  choice  of  future  and  present  packages,  careful  thought  must  be  given 
to  the  input  parameters  and  their  ranges. 

In  the  present  study,  the  analysis  parameter  included  the  solder  relative  shear 
displacement.  A  range  of  1.0  to  3.0  microns  is  chosen  for  the  relative  shear 
displacement  ( u ).  Due  to  the  size  of  the  package  and  FE  analysis  results,  the  solder 
joint  bending,  normally  represented  by  a  rotation  of  the  top  pad,  was  neglected.  In 
addition,  work  due  to  local  mismatch  at  each  interface  and  relative  axial 
displacement  is  ignored  and  set  equal  to  the  solder  free  expansion  since  it  is 
believed,  and  past  simulations  have  shown,  that  the  global  relative  shear  and  bending 
were  dominant. 

Solder  height  and  volume  parameters  were  chosen  for  the  study.  These 
parameters  are  normally  considered  design  parameters  along  with  solder  pad 
diameter,  however,  manufacturing  variations,  which  include  warpage  and  solder 
volume  variation,  can  affect  the  height  of  the  solder  joint.  Therefore,  the  present 
research  will  classify  these  parameters  as  design/manufacturing  parameters.  Three 
different  values  (levels)  of  the  variables  were  selected.  These  values  were  assumed 
to  be  uniformly  distributed  about  the  nominal  value  in  a  range  of  ±20  percent  for 
volume  and  +18  percent  for  height.  This  choice  of  the  range  was  motivated  by  a 
belief  that  even  within  a  reasonably  small  range  around  the  nominal  value,  an 
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optimal  solution  other  than  the  nominal  solution  must  exist.  For  this  study,  the 
solder  upper  and  lower  pad  diameters  are  held  constant. 

One  pure  manufacturing  parameter,  pad  tilt  of  individual  solder  joints  during 
reflow  due  to  warpage,  was  chosen  for  the  study.  The  warpage  represented  by  the 
angle  6warr  was  assumed  to  have  a  range  of  0.0°  to  3.0°.  The  parameters  (factors) 
and  levels  are  shown  in  Table  4. 1 . 


Table  4.1  Factors  and  levels. 


Factor 

Level  1 

Level  2 

Level  3 

A:  Volume  (V/V0) 

0.80 

B:  Height  (H/Ho) 

1.18 

C:  Pad  Tilt  due  to  Warpage  (deg) 

0.00 

IHSTOH 

D:  Relative  Shear  Displacement  (pm) 

1.00 

2.00 

3.00 

Vo  =  2.5  E8  pm3 

Ho  =  550  pm 

It  is  parenthetically  noted  that  to  build  a  global  approximation  in  a  problem 
with  a  large  number  of  variables,  a  partial  (fractional)  factorial  experiment  such  as 
an  orthogonal  array  of  Taguchi  (Peace,  1993)  may  be  more  appropriate  for  economy 
of  analysis.  A  designed  experiment  using  a  Taguchi  L27  orthogonal  array  based  on 
three  levels  of  design,  analysis,  and  manufacturing  variables,  as  shown  in  Table  4.2, 
was  used  to  capture  the  response  of  the  solder  joint. 
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Table  4.2  The  fractional  designed  experiment  used  to  build  the 
regression  model  for  solder  work. 


Factor 

A 

B 

C 

D 

Run 

Volume  Height  Warpage  (deg)  Shear  (mm) 

Work  (N-mm) 

1 

0.80 

0.82 

0 

1 

1.8380 

2 

0.80 

0.82 

1.5 

1 

1.8429 

3 

0.80 

0.82 

3 

1 

1.8565 

4 

0.80 

1.00 

0 

2 

4.1110 

5 

0.80 

1.00 

1.5 

2 

4.1180 

6 

0.80 

1.00 

3 

2 

4.1280 

7 

0.80 

1.18 

0 

3 

5.9770 

8 

0.80 

1.18 

1.5 

3 

5.9670 

9 

0.80 

1.18 

3 

3 

5.9360 

10 

1.00 

0.82 

0 

2 

5.4920 

11 

1.00 

0.82 

1.5 

2 

5.4590 

12 

1.00 

0.82 

3 

2 

5.5050 

13 

1.00 

1.00 

0 

3 

8.1220 

14 

1.00 

1.00 

1.5 

3 

8.0030 

15 

1.00 

1.00 

3 

3 

8.2080 

16 

1.00 

1.18 

0 

1 

1.2311 

17 

1.00 

1.18 

1.5 

1 

1.2392 

18 

1.00 

1.18 

3 

1 

1.2236 

19 

1.20 

0.82 

0 

3 

10.3420 

20 

1.20 

0.82 

1.5 

3 

10.5040 

21 

1.20 

0.82 

3 

3 

10.3260 

22 

1.20 

1.00 

0 

1 

1.6764 

23 

1.20 

1.00 

1.5 

1 

1.6887 

24 

1.20 

1.00 

3 

1 

1.6738 

25 

1.20 

1.18 

0 

2 

3.9680 

26 

1.20 

1.18 

1.5 

2 

3.9860 

27 

1.20 

1.18 

3 

2 

3.9890 

A  non-linear  elastic-plastic  analysis  was  performed  on  each  of  the  twenty- 
seven  settings  in  the  designed  experiment  using  a  new  Surface  Evolver  script  that 
produces  a  3-D  FE  surface  mesh  of  Evolver’ s  solder  profile  for  import  into  the 
IDEAS  finite  element  code.  The  IDEAS  code  is  used  to  generate  the  3-D  finite 
element  solid  mesh.  The  3-D  finite  element  mesh  is  analyzed  using  the  ABAQUS 
finite  element  code.  A  typical  finite  element  model  of  the  solder  joint  generated  by 
using  this  methodology  is  shown  in  Figure  4.1.  The  model  incorporates  over  1700 
tetrahedral  3-D  solid  elements.  The  accelerated  test  condition  consists  of  a  two 
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minutes  linear  temperature  loading  ramp  with  no  dwell  period  and  the  temperature 
range  is  from  0°C  to  100  °C.. 


Figure  4.1  A  typical  FE  model  of  the  solder  joint. 

A  linear  regression  model  (Equation  4.1)  was  built  using  the  commercial 
software  DOE  KISS  to  relate  the  work  quantities  with  the  four  input  parameters. 

W  =  4.75597  +  0.50607  *  A  -  0.83398  *  B  +  0.0049 1  *C  + 

3. 1 195 1  *  D  -  0.32930  *  A  *  B  +  0.5 1 520  *  A  *  D  -  0.36337  *  B  *  D 

The  fit  is  accurate  as  measured  by  an  R2  value  of  0.9984. 

The  relative  displacement  can  be  identified  for  each  solder  joint  after 
temperature  cycling.  This  can  be  used  to  calculate  the  stress/strain  distributions  and 
the  associated  inelastic  dissipation  terms  with  respect  to  temperature  loading  levels, 
rates,  and  frequencies.  For  this  study,  the  calculation  includes  only  elastic  and  rate- 
independent  plastic  behavior  of  solder.  The  output  is  used  for  fatigue  life  estimation. 

4.3.3  Solder  Joint  Fatigue  Life  Calculation 

Several  fatigue  indicators  have  been  used  to  predict  fatigue  failure  in  solder 


joints.  These  include: 
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a.  Maximum  plastic  strain  (Corbin,  1993) 

b.  Plastic  strain  range  (Nagaraj  and  Mahalingam,  1993) 

c.  Inelastic  strain  energy  density  (Darveaux,  et  al.,  1995) 

d.  Strain  energy  partitioning,  (Barker,  et  al.,  1993). 

Failure  of  ductile  solder  material  arises  from  repeated  nonlinear  deformation. 

Previous  studies  have  shown  that  strain  energy  density  gives  a  very  good  correlation 
to  solder  crack  growth  data  (Wong,  et  al,  1999).  An  empirically  derived  formula 
determined  by  Darveaux,  et  al.  (1995)  is  modified  and  used  in  the  current  study  to 
estimate  fatigue  life.  The  formula  relates  the  amount  of  inelastic  strain  energy  density 
developed  during  one  temperature  cycle  (AW)  to  the  number  of  cycles  needed  to 
induce  solder  failure,  Nf,  and  is  described  below. 


AC  oc 


AW 


4.2 


where  AW®  is  the  inelastic  strain  energy  density  at  the  solder/package  or 
solder/PWB  interfaces. 

To  determine  Nf  from  equation  4.2,  we  need  to  estimate  AW®.  For 
Accelerated  Temperature  Cycling  (ATC),  the  solder  inelastic  strain  energy  density  at 
the  end  of  the  cycle  and  the  inelastic  strain  energy  density  per  cycle  for  each 
element,  AWet,  can  be  obtained  from  the  FE  analysis  of  the  solder  model.  Since 
failure  usually  occurs  at  the  solder/package  or  solder/PWB  interface,  it  is  reasonable 
to  assume  the  inelastic  strain  energy  density  is  distributed  equally  between  the  two 
interfaces.  The  interfacial  strain  energy  density  per  cycle  is 


A^nt 


2VcAflr 


4.3 
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where: 

A Wei:  Inelastic  Strain  Energy  Density/Element 
Vei :  Element  Volume 

Vchar ■  Characteristic  Volume  at  the  solder/package  or  solder/PWB 

interface  contributing  to  flaw  growth 
Vchar  =  Pod  Area  *  Characteristic  Length 
In  this  study,  VChar  is  the  same  for  each  solder  joint,  therefore, 

— - *SF  = - - - *SF  4.4 

AWeIVel  Total  Plastic  Energy  Dissipation 

where  SF  is  a  Scale  Factor  used  specifically  for  this  study  to  generate  whole  number 
quantities  for  failure  cycles.  For  this  study,  SF  was  set  equal  to  100,000.  Since  we 
only  considered  elastic  and  rate-independent  plastic  solder  behavior,  the  applied 
loading  on  the  solder  joints  consisted  of  a  modified  ATC  specified  by  Darveaux.  et  al. 
(1995):  two  minutes  linear  temperature  loading/unloading  ramps  with  no  dwell  period 
and  the  temperature  range  is  from  0°C  to  100  °C. 

It  is  reasonable  to  expect  that  any  quantity  such  as  the  inelastic  dissipation 
calculated  as  a  result  of  finite  element  analysis  is  dependent  on  finite  element  mesh 
size.  Therefore,  a  mesh  density  study  was  performed  to  determine  the  appropriate 
mesh  density  for  the  solder  model.  The  element  edge  lengths  were  varied  for  a 
nominal  solder  joint  under  a  constant  relative  shear  over  the  above  temperature  cycle. 
The  characteristics  of  the  nominal  solder  joint  are: 

a.  Height  =  550  pm 

b.  Volume  =  2.5  E8  pm3 
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c.  Upper/Lower  Pad  Diameter  =  600  |im 

d.  Relative  Shear  Displacement  =  2.0  |im 

The  total  inelastic  energy  dissipations  and  errors  were  plotted  and  are  shown  in  Figure 
4.2. 


Average  Number  of  Elements 


261000  53500  10000  1700  250 


— Plastic  Energy 

— Percent  Error  Between 
Successive  Edge  Lengths 
— Percent  Error  as  Compared  to 
Edge  Length  9.8 
- Energy  Trendline 


Figure  4.2  Mesh  density  study. 

This  study  caused  an  element  edge  length  of  36.2  Jim  (10,000  elements) 
being  chosen  for  the  solder  joint  models.  The  resulting  plastic  dissipation  percent 
error  as  compared  to  element  edge  length  of  9.8  fim  (261,000  elements)  is  6.24 
percent,  which  is  reasonable  for  this  study. 

A  linear  regression  model  (Equation  4.5)  was  built  using  DOE  KISS  to  relate 
solder  fatigue  life  with  the  four  input  parameters,  volume,  height,  warpage  tilt,  and 
relative  shear  displacement,)  for  the  2x2  hypothetical  warped  model. 
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Nf  =  5253.9  -  216.295  *  A  +  566.348  *B-  24.538  *  C  - 

2796.2  *  Z)  +  2.782  *  A*  B*  C  -  802.0  *  A*  B*  D-  4.5 

15.121*  A*C*D2 -12.87  *B*C*D 

The  accuracy  of  the  fit  as  measured  by  its  R  value  was  0.9813.  The  fatigue 
life  expression  is  obtained  from  the  original  Taguchi  L27  designed  experiment  based 
on  three  levels  of  design,  analysis,  and  manufacturing  variables  listed  in  Table  4.1  and 
the  solder  joint  FE  model  described  in  Figure  4.1.  Results  of  the  experiment  is  shown 
in  Table  4.4 


Table  4.3  The  fractional  designed  experiment  used  to  build  the 
regression  model  for  solder  fatigue  life. 

Factor  ABC  D 

Run  Volume  Height  Warpage  (deg)  Shear  (mm)  Cycles  to  Failure 


1 

0.80 

0.82 

0 

1 

8658 

2 

0.80 

0.82 

1.5 

1 

8681 

3 

0.80 

0.82 

3 

1 

8518 

4 

0.80 

1.00 

0 

2 

5308 

5 

0.80 

1.00 

1.5 

2 

5342 

6 

0.80 

1.00 

3 

2 

5325 

7 

0.80 

1.18 

0 

3 

3971 

8 

0.80 

1.18 

1.5 

3 

3965 

9 

0.80 

1.18 

3 

3 

3935 

10 

1.00 

0.82 

0 

2 

4268 

11 

1.00 

0.82 

1.5 

2 

4202 

12 

1.00 

0.82 

3 

2 

4207 

13 

1.00 

1.00 

0 

3 

3056 

14 

1.00 

1.00 

1.5 

3 

3094 

15 

1.00 

1.00 

3 

3 

3045 

16 

1.00 

1.18 

0 

1 

8718 

17 

1.00 

1.18 

1.5 

1 

8688 

18 

1.00 

1.18 

3 

1 

8688 

19 

1.20 

0.82 

0 

3 

2509 

20 

1.20 

0.82 

1.5 

3 

2396 

21 

1.20 

0.82 

3 

3 

2430 

22 

1.20 

1.00 

0 

1 

8203 

23 

1.20 

1.00 

1.5 

1 

8065 

24 

1.20 

1.00 

3 

1 

8110 

25 

1.20 

1.18 

0 

2 

5230 

26 

1.20 

1.18 

1.5 

2 

5238 

27 

1.20 

1.18 

3 

2 

5222 
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In  general,  the  ability  of  linear  regression  models  to  fit  highly  nonlinear  data 
is  limited  by  the  nature  of  the  assumed  regression  model.  For  example,  a  quadratic 
model  would  provide  a  poor  fit  over  a  wide  range  to  a  data  that  has  an  exponential 
form.  Since  it  is  difficult  to  know  the  functional  form  of  the  data  a  priori,  several 
regression  models  would  have  to  be  tried  before  a  suitable  model  can  be  chosen. 
The  automated  regression  analysis  technique  incorporated  into  DOE  KISS  affords  us 
an  opportunity  to  examine  different  interactions  within  a  design  to  determine  the 
most  appropriate  model.  However,  Deshpande,  et  al.  (1997)  and  Subbarayan,  et  al. 
(1996)  showed  artificial  neural  network  (ANN)  models  would  provide  a  better 
approximation  than  a  linear  model  since  the  problem  surfaces  in  fatigue  analysis  are 
highly  nonlinear.  This  is  due  to  the  mathematically  proven  ability  of  ANN  models  to 
fit  arbitrarily  complex  functions.  The  use  of  ANN  models  is  beyond  the  scope  of  the 
present  study  and  a  linear  regression  model  will  be  used  to  get  a  general  idea  on  how 
the  fatigue  life  varies  over  the  range  of  the  factors. 

4.3.4  Validation  of  Methodology 

A  hypothetical  2x2  area-array  electronic  assembly  is  used  to  validate  the 
methodology.  The  material  properties  for  this  assembly  are  the  same  as  those  used  for 
the  3-D  electronic  packages  in  Chapter  2.  In  this  assembly,  solder  volume  variation  is 
assumed  responsible  for  the  variations  in  shapes  of  solder  joints.  The  solder  volume  is 
varied  linearly  from  the  lower  left  comer  of  the  assembly  to  the  upper  right  comer  as 
shown  in  Figure  4.3.  By  varying  the  volume  in  this  manner,  three  out  of  the  four 
solder  joints  had  different  shapes.  The  shape  prediction  program  is  then  used  to 
predict  the  final  equilibrium  configuration  of  the  assembly.  The  equilibrium 
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configuration,  shown  in  Figure  4.3  contained  a  package  tilt  of  1.77  degrees.  A  full 
three-dimensional  finite  element  model  of  the  warped  assembly  was  also  built  and  is 
shown  in  Figure  4.4.  A  side  view  of  the  warped  assembly  showing  the  package  tilt  is 
shown  in  Figure  4.5.  The  assembly  is  viewed  from  the  lower  right  comer  to  the  upper 
left  comer.  The  cross-section  of  a  nominal  “ideal”  assembly  is  shown  in  Figure  4.6. 


Vol:  0.250  mmA3  Vol:  0.200  mmA3 

Hght:  0.555  mm  Hght:  0.490  mm 

Tilt:  1.77  deg  Tilt:  1.77  deg 


Vol:  0.300  mmA3  Vol:  0.250  mmA3 
Hght:  0.620  mm  Hght:  0.555  mm 

Tilt:  1.77  deg  Tilt:  1.77  deg 


Figure  4.3  Equilibrium  configuration 
with  package  tilt. 


Figure  4.4  3-D  model  of  the  2x2  assembly. 


Figure  4.5  Side  view  of  warped  assembly. 
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MOLD 

DIE 

SUBSTRATE 


PCB 

D1=D2  =  0.600 
Solder  Volume  =  0250  mmA3 
Solder  Height  =  0.558  mm 
All  dimensions  mm 


Figure  4.6  Cross-section  of  nominal  assembly. 


The  full  3-D  FE  model  of  the  2x2  assembly  was  decomposed  into  linear  and 
nonlinear  sub-structures,  package,  PCB  and  solder.  Superelements  were  generated 
from  the  3-D  FE  models  of  the  linear  sub-structures,  package  and  PCB,  and  using  the 
work  regression  model  for  the  non-linear  sub-structure,  solder,  decomposition  was 
used  to  analyze  the  assembly. 

4.3.5  Results  and  Discussion 

Table  4.5  shows  the  results  of  the  analysis  using  the  decomposed  solution  as 
compared  to  the  finite  element  solution  of  the  warped  and  nominal  assemblies.  Figure 
4.7  shows  the  decomposition  and  FE  relative  shear  displacement  results  and  errors  for 
the  warped  assembly.  The  2x2  warped  area-array  model  has  an  average  displacement 
error  of  12.57%  across  all  four  solder  joints  while  the  overall  average  relative  shear 
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and  axial  displacement  errors  were  5.06%  and  5.62%  respectively.  The  2x2  nominal 
area-array  model  has  an  average  displacement  error  of  12.41%  across  all  four  solder 
joints  while  the  overall  average  relative  shear  and  axial  errors  were  5.51%  and  4.70% 
respectively.  The  high  average  displacement  errors  can  be  attributed  to  our  neglecting 
the  work  due  to  the  local  thermal  mismatch  at  the  solder  interfaces.  We  showed  in 
Chapter  2  that  by  including  this  analysis  variable  we  were  able  to  reduce  the  average 
error  for  the  225  I/O  PBGA  package  by  over  thirty  percent  from  8.0%  to  5.5%.  It  is 
reasonable  to  expect  the  same  error  reduction  for  the  warpage  study,  however,  as 
previously  noted,  solder  reliability  is  based  on  the  relative  displacements  and  our 
relative  displacement  errors  were  very  low.  Due  to  symmetry,  only  one-quarter  of  the 
nominal  package  was  analyzed  using  FE  analysis. 

Table  4.4  Results  of  decomposition  analysis  of  warped  assembly. 


Displacement  Error  (%) 
u  v  w 

Relative  Displacement 
Error  (%) 

shear  axial 

Warped  Assembly 

Component 

PWB 

17.84  18.08  13.18 

4.86  4.69  16.77 

Average  12.57 

5.06  5.62 

Nominal  Assembly 

Component 

PWB 

17.86  18.06  13.21 

4.56  4.34  16.46 

Average  12.41 

5.51  4.70 

Decomp/Error:  2.010  |i.m/6. 12% 
3D  FE  w/  Package  Tilt:  2.141  |im 
3D  FE  Nominal:  2.13  pm 


Decomp/Error:  1.918  |im/5. 33% 
3D  FE  w/  Package  Tilt :  2.026  |im 
3D  FE  Nominal:  2.13  |im 


Decomp/Error:  2.108  (im/ 1.78% 
3D  FE  w/  Package  Tilt :  2.146  fim 
3D  FE  Nominal:  2.13  (im 


Decomp/Error:  2.009  (im/5.86% 
3D  FE  w/  Package  Tilt :  2.134  [im 
3D  FE  Nominal:  2.13  |im 


Figure  4.7  Relative  shear  displacements  and  errors. 


A  reliability  study  was  completed  to  study  the  affect  of  warpage  on  solder 
joint  reliability.  The  study  compared  results  of  the  decomposition  and  FE  modeling. 
Critical  to  this  study  was  the  effect  of  the  neglecting  package  tilt  on  solder  joint 
reliability.  To  accomplish  this,  the  package  tilt  parameters  in  the  solder  shape 
program  were  turned  off  and  the  package  reanalyzed  to  find  the  equilibrium 
configuration.  Figure  4.8  depicts  the  equilibrium  configuration  without  package  tilt. 
Additionally,  relative  shear  displacements  for  FE  model  without  package  tilt  were 
assumed  the  same  as  those  for  the  FE  model  with  package  tilt.  Results  of  the 
reliability  study  are  shown  in  Figures  4.9  and  4.10.  The  percent  errors  in  Figure  4.9 
are  based  on  the  decomposition  and  FE  results  of  the  2x2  hypothetical  model  with 


package  tilt. 
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Vol:  0.250  mmA3 
Hght:  0.554  mm 
Tilt:  0.00  deg 


Vol:  0.300  mmA3 
Hght:  0.554  mm 
Tilt:  0.00  deg 


Vol:  0.200  mmA3 
Hght:  0.554  mm 
Tilt:  0.00  deg 


Vol:  0.250  mmA3 
Hght:  0.554  mm 
Tilt:  0.00  deg 


Figure  4.8  Equilibrium  configuration 
without  package  tilt. 


Decomp  w/  Package  Tilt:  5251/5.0% 
3D  FE  w/  Package  Tilt:  5000 
3D  FE  w/o  Package  Tilt:  5010 
3D  FE  Nominal:  5080 


Decomp  w/  Package  Tilt:5704/3.5% 
3D  FE  w/  Package  Tilt:  5510 
3D  FE  w/o  Package  Tilt:  5160 
3D  FE  Nominal:  5080 


Decomp  w/  Package  Tilt:  4774/0.01% 
3D  FE  w/  Package  Tilt:  4810 
3D  FE  w/o  Package  Tilt:  5290 
3D  FE  Nominal:  5080 


Decomp  w/  Package  Tilt:  5251/4.6% 
3D  FE  w/  Package  Tilt:  5020 
3D  FE  w/o  Package  Tilt:  5040 
3D  FE  Nominal:  5080 


Figure  4.9  Results  of  reliability  study  (cycles  to  failure). 
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Reliability  Study 


— Decomposed  w/  tilt  — ■—  FE  w/  tilt  -±—  FE  w/o  tilt  — Nominal 
Figure  4.10  Cycles  to  failure  plot. 

As  seen  in  Figure  4.10,  for  this  hypothetical  package,  failure  to  include 
package  tilt  in  determining  the  equilibrium  of  warped  electronic  assemblies  can  result 
in  errors  in  reliability  estimates  by  as  much  as  10  percent.  The  reason  for  this  error  is 
volume  variation  in  the  assembly  caused  the  package  to  tilt,  which  directly  affected 
the  standoff  height  of  the  solder  joints.  Chan,  et  al.  (1997)  showed  in  their  study  that 
fatigue  life  is  strongly  dependent  on  the  height  of  solder  joint.  This  is  a  critical 
finding  because  as  electronic  packages  continue  to  grow  in  size,  warpage  and  its  effect 
of  package  tilt  will  increase  in  importance,  and  ignoring  tilt  would  lead  to  large  errors 
in  reliability.  For  example,  an  1140  I/O  Flip  Chip  Ball  Grid  Array  (FCBGA) 
electronic  package  shown  in  Figure  4.11  has  a  diagonal  length  of  26.9  mm.  The 
package  has  a  500  micron  nominal  solder  joint  height  and  the  solder  ball  array  is 
peripherally  distributed.  A  package  of  this  type  would  only  need  a  1.065  degree  tilt  at 
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the  center  to  make  the  edge  of  the  package  incident  on  the  PCB.  We  must  note  the 
volume  variation  in  our  hypothetical  package  was  deliberately  chosen  to  be  extreme  to 
properly  demonstrate  the  methodology.  We  would  not  expect  to  find  such  a  volume 
distribution  in  “real-world”  electronic  packages,  however,  arbitrary,  non-symmetric 
warpage  of  PCBs  and  packages  is  a  major  concern  in  the  electronic  packaging 
industry,  and  we  showed  in  Chapter  3  warpage  of  this  type  could  cause  the  package  to 
tilt.  The  FCBGA  package  above  would  not  require  such  an  extreme  volume  variation 


or  arbitrary  warpage  to  achieve  a  1.065  degree  package  tilt. 


Figure  4.11  1 140  I/O  Flip  Chip  Ball  Grid  Array  Package. 

Lastly,  using  the  regression  model  for  solder  work,  Equation  4.1,  plots  of  the 
interactions  between  the  design/manufacturing  variables,  volume,  height  and  warpage, 
and  analysis  variable,  relative  shear  (Figures  4.12  -  4.14),  show  the  interactions 
between  height  and  shear  and  warpage  and  shear  to  be  negligible.  The  absence  of 
interaction  is  evident  in  the  parallel  running  of  the  lines.  There  is  a  small  interaction 
between  volume  and  shear,  however,  for  this  study  the  solder  volume  distribution  was 
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extreme  to  demonstrate  the  methodology.  As  stated  previously,  we  would  not  find 
such  a  volume  distribution  in  a  “real-work”  electronic  package.  This  suggests  that  the 
assembly  design/manufacturing  parameters  (solder  volume,  height,  and  warpage)  do 
not  significantly  influence  the  relative  displacements  on  the  solder  joints.  Further,  a 
comparison  of  the  FE  relative  shear  displacements  for  the  hypothetical  warped  and 
nominal  area-array  assemblies,  as  shown  in  Figure  4.7,  produced  a  maximum  relative 
shear  displacement  error  of  5.13  percent  at  solder  joint  #3.  That  is,  the  package  tilt 
does  not  significantly  alter  the  displacements  in  this  particular  package.  Therefore,  it 
is  reasonable  to  assume  that  the  analysis  of  nominal  area-array  assemblies  is  all  that  is 
needed  to  determine  the  assembly  equilibrium  configuration  during  temperature 
cycling.  The  relative  shear  displacements  from  these  analyses  are  sufficient  for 
subsequent  solder  joint  reliability  studies. 


+ 

0002  0.00021  0.00022  0.00023  0.00024  0.00025  0.00026  0.00027  0.00028  0.00029 

Volume  (cm3) 


0.0003 


Figure  4.12.  Interaction  plot  of  volume  vs.  shear. 
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Interaction  Plot  of  Warpage  vs  Shear 
Constants:  Volume  =  0.00025  cm3  Height  =  550  microns 
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Figure  4.14  Interaction  plot  of  warpage  vs.  shear. 

3-D  surface  plots  of  the  linear  regression  model  for  solder  fatigue  life 
(Equation  4.5)  provide  important  insight  into  the  effect  of  warpage  on  solder  joint 
reliability.  A  pad  tilt  caused  by  warpage  can  affect  reliability  in  two  ways:  due  to  a 
changed  standoff  height  (and  its  subsequent  effect  on  stresses)  or  due  to  the  change 
in  the  nature  of  singularity  at  the  joint/PCB  or  the  joint/component  interface.  A  plot 
of  solder  height  versus  pad  tilt  due  to  warpage  (Figure  4.15)  shows  solder  pad  tilt 
caused  by  warpage  independent  of  solder  standoff  height  does  not  affect  reliability. 
That  is,  the  nature  of  the  singularity  is  not  significantly  altered  by  the  pad  tilt.  A  plot 
of  solder  volume  versus  height  (Figure  4.16)  shows  that  as  the  shape  of  the  solder 
changes  from  convex  (truncated  sphere)  to  concave  (hour-glass)  the  fatigue  life  of 
the  solder  joint  increases.  The  main  reason  for  this  change  in  reliability  versus 
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solder  shape  change  is  because  the  inelastic  dissipation  is  concentrated  at  the 
location  having  the  smallest  cross-section.  In  addition,  convex  (truncated  sphere) 
solder  joint  have  higher  strain  energy  density  due  to  the  acute  contact  angle  at  the 
solder/package  and  solder/PCB  interfaces  causing  (elastically)  singular  stress  fields. 
As  shape  changes  from  convex  to  concave  (hour-glass  shaped)  the  interfacial  contact 
angle  increases  reducing  the  stress  concentrations  (Deshpande,  et  al,  1997).  This  can 
be  seen  in  Figures  4.17  -  4.20  which  show  the  plastic  dissipation  contours  of  a 
solder  joint  with  a  prescribed  volume  (2.5  E-4  cm3),  fixed  relative  shear 
displacement  (2.0  microns),  and  prescribed  pad  diameters  (600  microns).  The  height 
of  the  solder  joint  varies  from  500  to  1100  microns,  as  is  necessary  to  produce 
convex  and  concave  shapes.  For  a  convex  (truncated  sphere)  shaped  solder  joint 
(Figure  4.17)  with  low  interfacial  contact  angle  (high  curvature),  the  highest  plastic 
dissipation  value  is  approximately  0.3661  MPa  and  is  located  at  the  pad  interfaces. 
For  a  barrel-shaped  solder  (Figure  4.18),  with  increased  interfacial  contact  angle,  the 
highest  strain  energy  density  value  is  approximately  0.3457  MPa  and  is  located  at 
the  pad  interfaces.  For  a  cylindrical  shaped  solder  (Figure  4.19),  with  a  further 
increase  in  interfacial  contact  angle,  the  highest  strain  energy  density  value  is 
approximately  0.1740  MPa  and  is  located  at  the  pad  interfaces.  Lastly,  for  a  concave 
(hour-glass)  shaped  solder  (Figure  4.20),  with  maximum  interfacial  contact  angle 
(low  curvature),  the  highest  strain  energy  density  value  is  approximately  0.1153 
MPa  and  is  located  close  to  the  center  of  the  solder  joint  where  the  cross-section  area 
is  smallest.  This  study  confirms  the  results  of  the  studies  by  Ju,  et  al.  (1994,  Chan,  et 
al.  (1997),  and  Deshpande,  et  al.  (1997). 
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Cycles  to  Failure 

Surface  Plot  of  Height  vs  Warpage 
Constants:  Volume  =  0.00025  cm3  Shear  =  2  microns 


Warpage  Tilt  (deg) 

□  4000-5000  E  5000-6000 


Figure  4.15  Cycles  to  failure  surface  plot  of  height  vs.  warpage. 


Barrel-Shaped  Solder 
Joint 


| Concave  Solder  Joint] 


Cycles  to  Failure 
Surface  Plot  of  Volume  vs  Height 
Constants:  Warpage  =  0  deg  Shear  =  2  microns 


|  □  4000-5000  H  5000-6000  ■6000-7000  | 


Figure  4.16  Cycles  to  failure  surface  plot  of  volume  vs.  height. 
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For  the  2x2  warped  assembly,  the  shape  of  the  solder  joint  ranged  from 
concave  to  barrel-shaped  (See  Figure  4.16).  Lastly,  Figure  4.21  shows  how  a 
combined  bowl-type  warpage  of  the  assembly  can  be  enhanced.  The  dotted  line 
indicates  the  cycles  to  failure  path  for  an  arch-type  warpage.  The  fatigue  life  ranges 
from  8000-10000  cycles  for  solder  joints  located  close  to  the  neutral  point  (low 
relative  shear)  of  the  assembly  to  2000-4000  cycles  for  solder  joints  located  far  away 
from  the  assembly  neutral  point  (high  relative  shear).  The  solid  line  indicates  the 
cycles  to  failure  path  for  a  combined  bowl-type  assembly  warpage.  The  fatigue  life 
ranges  from  6000-8000  cycles  for  solder  joints  located  close  to  the  assembly  neutral 
point  to  0-2000  cycles  for  joints  located  far  away  from  the  neutral  point.  Therefore, 
depending  on  the  type  of  package  warpage,  the  fatigue  life  at  solder  joints  with  high 
relative  shear  can  change  by  over  200  percent.  This  agrees  with  the  findings  of 
Chan,  et  al.  (1997). 
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Figure  4.21  Cycles  to  failure  surface  plot  of  height  vs.  shear. 
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4.4  Conclusions 

A  novel  procedure  is  developed  for  analysis  of  area-array  solder  joint 
reliability.  The  model  consists  of  three  steps:  1)  calculation  of  solder  shapes  using 
novel  package  configuration  prediction  methodology,  2)  calculation  of 
stresses/strains  and  cyclic  inelastic  dissipations  using  a  decomposition  modeling 
technique,  and  3)  estimation  of  solder  joint  fatigue  life  using  a  regression  model  for 
fatigue  life.  Fatigue  life  calculation  was  based  on  a  modified  Darveaux  relationship 
(Darveaux,  et  al.,  1995).  It  is  shown  that  the  decomposed  solution  methodology 
results  in  an  overall  accuracy  loss  in  relative  shear  displacements  of  approximately 
5.06  percent  when  compared  with  a  full  finite  element  model.  This  displacement  error 
leads  to  an  error  in  fatigue  life  of  approximately  3.28  percent  due  to  incorrectly 
calculated  inelastic  dissipation.  A  comparison  of  the  reliability  of  a  package  with 
warped  and  non-warped  joints  reveals  package  tilt  due  to  solder  volume  variation 
and/or  non-symmetric  warpage  can  change  the  stand-off  height  of  solder  joints 
significantly  affecting  solder  joint  reliability.  However,  we  found  that  individual 
solder  pad  tilts,  independent  of  their  effect  on  standoff  height,  do  not  affect  the 
reliability  of  solder  joints. 


Chapter  5 


Conclusions  and  Recommendations 

5.1  Conclusions 

During  the  course  of  the  research  presented  in  this  thesis,  a  novel 
methodology  for  assessing  the  reliability  of  warped  solder  joints  in  area-array 
electronic  packages  was  developed.  The  methodology  is  incorporated  into  the  three 
basic  and  critical  steps  to  determine  package  reliability: 

1 .  calculate  solder  j  oint  shapes , 

2.  calculate  thermal  stress/strain  distribution  in  the  most  susceptible  joint 
after  solidification  and  during  thermal  cycling,  and 

3.  assess  solder  joint  reliability. 

For  step  one,  a  novel  methodology  was  developed  to  determine  the  final 
equilibrium  configuration  of  an  area-array  package.  The  novel  techniques  in  this 
methodology  include; 

=>  First  application  of  homogeneous  transformation  theory  to  predict  the 
equilibrium  configuration  of  electronic  packages 
♦  First  use  of  rigid  rotations  of  package 
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=>  Methodology  is  capable  of  predicting  equilibrium  in  presence  of  arbitrary 
warpage 

For  step  two,  a  decomposed  analysis  procedure  is  successfully  implemented 
to  determine  the  deformation  of  area-array  assemblies  with  warped  solder  joints,  a 
novel  application  of  this  procedure.  In  addition,  analysis  of  2-D  and  3-D 
hypothetical  and  “real-world”  warped  and  non-warped  electronic  packages  produced 
the  following: 

=>  Refined  model  for  solder  response  through  additional  basis  terms 
=>  Successfully  demonstrated  methodology  on  “real-world”  electronic 
package  with  non-warped  solder  joints 
=>  Successfully  demonstrated  methodology  on  hypothetical  electronic 
package  with  warped  solder  joints 
=>  Provided  insight  into  decomposed  approach 

♦  Displacement  versus  relative  displacement  errors 
=>  Provided  insight  into  energy  contributions  to  relative  displacement  errors 
=>  Led  to  assumption  that  the  analysis  of  a  nominal  area-array  assembly  is 
all  that  is  required  to  obtain  the  displacements  required  to  assess 
reliability  of  individual,  warped  solder  joints. 

For  step  three,  our  study  showed  the  following: 

=>  Solder  pad  tilt  in  a  range  of  0  to  3  degrees,  independent  of  its  effect  on 
solder  stand-off  height,  has  a  negligible  effect  on  solder  joint  reliability 
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=>  Package  tilt  due  to  solder  volume  variation  and/or  arbitrary  package  or 
PCB  warpage  can  change  the  stand-off  height  of  solder  significantly 
affecting  solder  joint  reliability 

♦  We  found  a  reliability  calculation  error  as  high  as  10  percent  for  the 
2x2  hypothetical  array 

5.2  Recommendations 

The  results  of  the  research  successfully  demonstrate  the  power  of  the 
decomposition  analysis  procedure  and  the  novel  methodology  for  determining  the 
equilibrium  configuration  of  area-array  packages  during  solder  reflow.  As  with  any 
new  technique,  additional  research  must  be  accomplished  before  the  technique  can 
be  accepted  as  state-of-the-art  and  become  an  industry  standard  analysis  technique 
like  finite  element  analysis.  The  following  is  a  list  of  recommendations  for  future 
studies  related  to  the  work  presented  in  this  thesis: 

5.2.1  Decomposition 

=>  Generate  solution  for  high  relative  displacement  errors  attributed  to 
superelement  energy  contribution 

=>  Incorporate  time-dependent  (creep)  behavior  into  solder  work  model 
=>  Demonstrate  methodology  on  a  “real-world”  warped  electronic  package 

♦  Verify  assumption  that  the  decomposed  analysis  of  nominal  area- 
array  assemblies  is  all  that  is  needed  to  determine  package 
equilibrium  during  temperature  cycling 
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-  This  is  based  on  the  low  relative  displacement  errors  between  FE 
analysis  of  hypothetical  assemblies  with  warped  and  non-warped 
solder  joints  and  the  negligible  manufacturing  and  analysis 
variable  interactions  discovered  during  regression  analysis  of 
designed  experiments 

=>  Use  decomposition  to  determine  optimal  construction  of  area-array 
packages 

=>  Produce  Windows/Unix  based  user-friendly  program  utilizing  this 
methodology 

5.2.2  Solder  Shape  Prediction 

=>  Incorporate  pad  off-set  into  shape  prediction  modeling 

=>  Demonstrate  procedure  on  “real-world”  electronic  package 

=>  Use  solder  shape  prediction  program  as  a  tool  to  predict  solder  joint  yield 
during  solder  reflow 

5.2.3  Fatigue  Life  Calculation 


=>  Assess  reliability  of  warped  solder  joints  in  the  presence  of  pad-offset 
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